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Abstract 

The  key  to  Direct  Simulation  Monte  Carlo  (DSMC)  is  decoupling  of  particle 
motion  and  particle  collisions.  Particles  within  each  cell  are  randomly  chosen  as  col¬ 
lision  partners,  and  the  collision  is  then  accepted  or  rejected  by  comparing  collision 
criteria  to  a  random  number.  In  the  Smoothed  Accept/Reject  (SAR)  algorithm,  the 
accept/reject  criteria  is  altered:  rather  than  a  binary  function  of  rejection  or  accep¬ 
tance,  collisions  can  be  partially  accepted  with  a  linear  weighting  between  zero  and 
one.  The  partial  acceptance  is  based  on  a  band  around  the  original  accept/reject 
criteria  defined  as  a  percentage  of  the  collision  criteria,  which  is  called  e.  The  weight¬ 
ing  is  used  in  sampling  the  particles  in  order  to  calculate  the  macroscopic  flowfield 
parameters.  Previous  work  included  comparisons  to  experimental  data  using  inverse 
shock  thickness,  the  results  of  which  implied  a  relationship  for  the  appropriate  value 
of  e  as  a  function  of  Mach  number,  and  this  is  explored  in  the  present  work.  Addi¬ 
tionally,  2-dimensional  experimental  data  is  used  for  comparison  between  SAR  and 
the  original  algorithm.  Velocity  distributions  of  the  particles  are  examined  for  all 
algorithms  and  compared  to  experimental  data  to  determine  the  effect  of  the  SAR 
algorithm  at  a  microscopic  level.  Both  of  the  1-dimensional  comparisons  to  experi¬ 
ment  shows  that  SAR  provides  results  that  best  matches  the  experimental  data.  All 
of  the  comparisons  to  experiment  show  a  Mach  dependency  that  has  previously  been 
noted,  and  the  dependency  was  defined  for  the  normal  shock  simulations.  DSMC 
does  adequately  simulate  the  nonequilibrium  within  the  cells  at  a  high  Mach  number 
through  the  shock,  but  SAR  does.  The  SAR  algorithm  models  the  flowfield  in  the 
shock  better  than  DSMC  through  a  change  in  the  collision  rate  and  particle  sampling 
methods,  which  allows  for  a  more  accurate  simulation. 
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Improved  Collision  Modeling  for  Direct  Simulation 

Monte  Carlo  Methods 

I.  Introduction 

1 . 1  Motivation 

The  United  States  Air  Force  (USAF)  has  twelve  core  functions  that  are  inte¬ 
gral  to  its  role  in  the  defense  of  our  nation.  Of  these  twelve  functions,  hypersonics 
plays  a  direct  role  in  five:  nuclear  deterrence  operations,  air  superiority,  space  supe¬ 
riority,  global  precision  attack,  and  rapid  global  mobility.  The  ability  to  understand 
hypersonic  flows  and  apply  the  research  to  the  design  of  hypersonic  vehicles  increases 
access  to  space  and  allows  a  greater  range  and  response  time  of  air  and  space  vehi¬ 
cles  throughout  the  world.  Recently,  the  new  report  on  Technology  Horizons  by  Dr. 
Werner  Dahm,  the  Chief  Scientist  of  the  Air  Force  [2],  emphasized  important  future 
capabilities  that  depend  on  hypersonics  research,  including  prompt  strike  systems 
that  would  cruise  at  Mach  6  and  a  two-stage-to-orbit  reusable  space  launch  vehicle, 
both  of  which  will  fly  through  a  rarefied  atmosphere.  Therefore,  it  is  imperative  for 
the  Air  Force  to  accurately  model  and  predict  the  behavior  of  these  vehicles  within 
the  rarefied  hypersonic  regime. 

At  higher  altitudes,  starting  around  100  km  [3],  the  air  is  less  dense  and  there¬ 
fore  the  Knudsen  number  (Kn)  is  higher.  The  Knudsen  number  is  a  measure  of  the 
rarefaction  of  a  gas,  and  will  be  explained  more  thoroughly  in  the  following  chap¬ 
ter.  Continuum  solvers  that  use  Navier-Stokes  or  Euler  equations  cannot  accurately 
predict  the  gas  behavior  in  high  Knudsen  number  regimes.  Figure  1.1  illustrates  the 
Knudsen  number  limits  for  the  continuum  solvers. 

Euler  based  solvers  are  used  in  the  inviscid  limit  (when  viscosity  is  at  or  very 
near  zero)  while  Navier-Stokes  based  solvers  can  be  used  for  a  simulation  with  a 
Knudsen  number  of  less  than  0.1.  The  continuum  assumption  breaks  down  at  about 
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Figure  1.1:  Knudsen  Number  Limits  on  Mathematical  Models  [4] 

Kn=0.1,  and  the  flow  can  be  considered  rarefied  above  that  value.  Direct  Simulation 
Monte  Carlo  (DSMC)  is  a  stochastic  method  which  utilizes  the  Monte  Carlo  statistical 
model  to  simulate  gas  behavior,  which  is  very  useful  for  these  rarefied  atmosphere 
hypersonic  simulations,  and  has  been  used  for  this  purpose  for  decades.  Notably, 
DSMC  has  been  used  for  simulations  of  the  space  shuttle  reentry  [3,5].  DSMC  has 
also  compared  well  to  many  wind  tunnel  tests,  which  have  been  used  to  enhance 
understanding  of  hypersonic  flows  [6,7].  Additionally,  DSMC  has  been  heavily  used 
in  the  study  of  microflows  [7].  DSMC  will  be  explained  in  more  detail  in  Chapter  2. 

There  are  two  other  methods  which  are  used  to  predict  behavior  in  a  rarefied 
condition.  The  first  is  molecular  dynamics  (MD),  which  is  a  direct  simulation  of 
each  particle,  and  therefore  has  generally  been  limited  to  nanoscale  simulations  [8]. 
MD  is  more  computationally  expensive  than  DSMC,  but  there  has  been  research  to 
improve  the  computational  efficiency  of  MD  in  the  last  couple  of  years  so  MD  can 
be  used  more  in  dilute  gas  simulations  [8].  MD  is  a  deterministic  simulation  of  the 
system  whereas  DSMC  is  a  statistical  model  of  the  system,  which  requires  simulated 
particles  representing  a  large  number  of  real  particles  in  order  to  get  accurate  results. 
The  methodology  of  DSMC  makes  it  more  computationally  efficient  compared  to 
MD,  even  though  results  of  DSMC  are  very  similar  to  MD.  The  other  method  is 
direct  numerical  simulation,  generally  referred  as  a  Boltzmann  solver.  Boltzmann 
solvers  are  also  computationally  expensive,  but  they  can  be  coupled  with  continuum 
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solvers  to  reduce  computational  cost  in  cells  that  have  a  low  local  Knudsen  number  [9] . 
Therefore  DSMC  is  a  just  as  accurate,  but  more  cost  effective  method  for  hypersonic 
method  simulations. 

There  are,  however,  some  problems  with  the  original  DSMC  algorithm  within 
areas  of  nonequilibrium  such  as  shocks  and  boundary  layers.  The  collision  rate  is  not 
high  enough  to  allow  for  proper  equilibration  compared  to  experimental  data.  A  few 
collision  models,  one  of  which  is  discussed  in  more  detail  in  Chapter  2,  were  created 
in  order  to  better  simualte  the  flow  in  nonequilibrium  conditions.  These  methods 
change  the  collision  rate,  and  also  alter  the  transport  properties,  which  only  occur  in 
nonequilibrium.  The  simple  code  that  this  project  uses  has  one  method  to  correct  the 
problem,  but  there  are  more  methods  that  have  been  developed  over  the  years.  The 
scope  of  this  project  is  not  to  evaluate  the  newer  collision  models,  but  to  investigate 
a  new  accept/reject  method  that  may  correct  DSMC  in  another  manner. 

Over  the  past  3  years,  a  modified  DSMC  algorithm  known  as  Smooth  Accep¬ 
t/Reject  (SAR)  was  developed  by  Greendyke  et  al  [10]  that  may  give  more  accurate 
results.  The  previous  research  investigated  a  potential  method  that  would  reduce  the 
convergence  time  for  DSMC  simulations.  The  research  showed  that  convergence  was 
not  effected  by  SAR,  but  the  flowfield  results  were  different.  The  previous  work  also 
showed  that  DSMC  is  not  always  accurate  at  higher  Mach  numbers.  A  one  dimen¬ 
sional  DSMC  shock  program  predicts  a  thicker  normal  shock  wave  when  comparing 
to  experimental  data  [10, 11]  for  velocities  above  Mach  3.  Further  discussion  of  the 
previous  work  can  be  found  in  the  next  chapter. 

The  purpose  of  this  thesis  is  to  further  refine  the  technique  used  in  the  SAR 
modification  and  to  better  understand  why  SAR  produces  the  results  that  have  been 
seen  thus  far.  The  new  technique,  once  understood  and  applied,  can  be  used  for 
development  of  aircraft  that  operate  within  the  rarefied  regime  at  hypersonic  veloc¬ 
ities.  The  methodology  of  the  project  includes  comparison  of  the  new  algorithm’s 
results  to  experiment  and  theoretically  understanding  the  effects  of  the  algorithm. 
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Additionally,  the  velocity  distribution  of  specified  cells  will  be  compared  to  the  speed 
distribution,  as  derived  from  the  Maxwellian  Distribution.  Mach  dependency  will  also 
be  investigated,  as  a  proportional  relationship  has  been  noted  in  previous  studies  [11]. 
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II.  Background 

In  order  to  have  an  understanding  of  the  project,  one  must  first  understand  the 
derivation  and  application  of  the  Boltzmann  equation.  The  Boltzmann  equation  can 
be  used  to  calculate  gas  flow  in  any  condition,  incuding  in  rarefied  and  nonequilibrium 
conditions.  Rarefaction  will  be  further  discussed  in  Section  2.2.  Nonequilibrium  can 
be  defined  as  a  perturbation  from  an  equilibrium  condition,  and  when  the  departure 
from  equilibrium  is  large,  the  conservation  equations  which  were  derived  using  a 
small  perturbation  from  equilibrium  are  incomplete.  The  moments  of  the  Boltzmann 
equation,  which  are  found  by  multiplying  the  Boltzmann  equation  by  a  function 
of  velocity  and  integrating  over  velocity  space,  must  be  applied  in  order  to  fully 
describe  the  fluid  flow  in  nonequilibrium  [12].  The  Boltzmann  equation  must  then 
be  simulated  in  order  to  calculate  the  flowfield,  which  can  be  done  using  DSMC. 
DSMC  uses  statistical  modeling  to  predict  the  collisional  behavior  of  a  gas  using  a 
Monte  Carlo  scheme  and  then  calculating  the  expected  motion  through  the  use  gas 
kinetics.  DSMC  requires  knowledge  of  kinetic  theory,  therefore  the  next  two  sections 
will  include  a  description  of  both  the  Boltzmann  equation  and  the  associated  kinetic 
theory  before  DSMC  is  explained. 

2.1  The  Boltzmann  Equation  and  Statistical  Mechanics 

The  Boltzmann  equation  “describes  the  rate  of  change,  with  respect  to  position 
and  time,  of  the  distribution  function.”  [12]  The  distribution  function  being  discussed 
is  the  velocity  distribution  function  of  the  particles  within  the  system.  In  otherwords, 
the  Boltzmann  equation  describes  the  molecular  motion  of  a  system,  which  can  be  used 
to  determine  the  overall  behavior  of  that  system.  Molecular  motion  can  be  described 
using  velocity  space;  particles  can  occupy  velocity  space  just  as  they  occupy  physical 
space.  Just  as  a  particle  has  a  position  in  physical  space  that  can  be  described  with 
coordinates  in  a  reference  frame,  a  particle  in  velocity  space  can  be  described  in  a 
coordinate  system.  The  velocity  element  that  the  particles  reside  in  can  be  defined 
as:  Ci  +  dci ,  c2  +  dc2,c3  +  dc3,  where  ci,c2,  and  c3  represent  the  three  dimensional 
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components  of  the  particle’s  velocity.  The  number  of  particles  within  the  element  can 
be  found  by: 

dN  =  N  f  {c)dcidc2dc%  (2.1) 

where  N  refers  to  the  total  number  of  particles,  and  f(c)  refers  to  the  velocity  func¬ 
tion.  Just  as  a  spatial  volume  can  be  defined  by  multiplying  the  three  dimensions  of 
elements  that  make  up  the  volume,  dVx  =  dxidx2dx3,  the  volume  in  velocity  space 
is:  dVc  =  dc\dc2dc3  [12].  The  Boltzmann  equation  can  be  derived  by  examining 
molecules  of  a  particular  velocity  class,  q,  within  a  volume  in  physical  space  and  a 
volume  within  velocity  space.  A  velocity  class  refers  to  a  group  of  particles  occupy¬ 
ing  the  same  velocity  element.  The  Boltzmann  Equation  accounts  for  the  complete 
behavior  of  the  particle  within  the  velocity  class:  [12]: 

S[nf(cj)]  S[nf(cj)]  8[Fjnf(cj)\ 

St  Cj  Sxj  Scj 

Where  St  is  the  change  in  time,  n  is  the  number  density,  Sxj  is  the  change  in  position, 
Scj  is  change  in  velocity,  and  /(q)  is  the  velocity  function.  The  first  term  in  Equation 
2.2  is  the  convection  term,  the  second  term  is  the  flux  of  molecules  across  the  surfaces 
of  the  volume,  the  third  term  represents  flux  of  molecules  into  the  velocity  volume 
due  to  external  forces  such  as  gravity,  and  the  right  hand  side  of  the  equation  is 
the  collision  integral,  which  represents  the  rate  of  change  in  the  number  of  molecules 
in  the  velocity  class  due  to  collisions.  The  convection  term  simply  refers  to  the 
movement  of  the  particles  through  the  system.  The  flux  of  molecules  accounts  for  an 
open  system  where  particles  are  moving  in  and  out  of  the  domain.  The  right  hand 
side  of  the  equation,  known  as  the  collision  integral,  is  what  makes  the  Boltzmann 
equation  so  difficult  to  solve  analytically.  The  computations  are  very  expensive,  and 
only  very  simple  geometries  with  a  fairly  small  domain  have  been  solved  using  a  direct 
Boltzmann  solver.  The  collision  integral  requires  knowledge  of  the  velocity  states  of 
the  particles  before  and  after  collision  in  order  to  be  calculated  [13],  which  can  be 
done  using  kinetic  theory  and  statistical  mechanics. 


<w(q)] 
St 


(2.2) 
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The  collision  integral  describes  the  rate  of  collisions  between  two  particles  be¬ 
longing  to  different  velocity  classes  at  a  certain  point  on  the  sphere  of  influence  as 
referenced  from  the  velocity  vector  and  the  line  of  centers.  The  sphere  of  influence 
refers  to  the  spherical  distance  at  which  a  particle  will  begin  to  effect  another  parti¬ 
cle.  If  intermolecular  forces  are  being  neglected,  the  sphere  of  influence  is  the  radius 
of  the  particle,  but  the  sphere  of  influence  will  be  larger  if  intermolecular  forces  are 
being  considered.  For  the  purpose  of  this  discussion,  sphere  of  influence  could  refer 
generically  to  either  of  these  definitions.  Line  of  centers  refers  to  a  line  that  can 
be  drawn  by  connecting  the  center  of  the  two  colliding  particles.  The  two  colliding 
particles  belong  to  different  velocity  classes,  which  are  called  c  and  £.  These  particles 
will  collide  and  change  velocity  classes  due  to  the  collision,  while  other  particles  will 
have  collisions  that  cause  them  to  enter  into  the  two  velocity  classes.  These  two  types 
of  collisions  are  referred  to  as  depleting  and  replenishing  collisions,  respectively.  The 
collision  term  can  now  be  expressed  as  the  integral  of  the  velocity  functions  over  the 
volume  of  the  sphere  of  influence  of  the  £  particles: 


S[nf(cj)] 
St 


OO  27T  71-/2 


J  coll 


™2[/(c')/(0  -  f(ci)f(Ci)]gd2sin'ipcos'ipd'ipdedV<:  (2.3) 


— oo  0  0 


where  c'  and  (■  replenish  the  velocities  while  ct  and  Q  deplete  the  velocity  classes,  g 
is  the  relative  velocity  between  the  particles,  d  is  the  radius  of  the  sphere  of  influence 
of  the  particle,  n  is  the  number  density,  and  e  and  ^  represent  angles  that  define 
the  location  of  the  collision  on  the  sphere  of  influence.  The  principle  of  reciprocity 
is  an  assumption  that  for  every  collision  that  depletes  from  a  velocity  class,  there  is 
a  collision  that  replenishes  from  another  velocity  class.  The  principle  holds  true  for 
equilibrium  cases  only.  Therefore,  in  equilibrium  the  collision  integral  equals  zero, 
which  means  the  integral  only  matters  for  nonequilibrium  conditions,  such  as  shocks 
or  boundary  layers  where  gradients  are  large  [12]. 

As  stated  before,  nonequilibrium  can  be  expressed  as  a  perturbation  away  from 
equilibrium,  and  occurs  in  two  areas  of  the  flowfield:  shock  layers  and  boundary  layers, 
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both  of  which  happen  in  hypersonics.  In  equilibrium,  particles  in  a  system  can  be 
distributed  into  energy  states  according  to  the  Boltzmann  Distribution  Function  [12]: 


NCje^r 

'LjCje^ 


(2-4) 


where  N  is  the  number  of  particles,  N*  is  the  number  of  j  macrostates,  Cj  is  the  num¬ 
ber  of  increments,  or  degeneracy,  k  is  the  Boltzmann  constant,  T  is  the  temperature, 
and  ej  is  the  energy.  The  Boltzmann  Distribution  Function  is  the  ratio  of  particles 
at  a  given  energy  level  over  the  sum  of  possible  energy  levels.  At  any  given  time,  the 
particles  may  not  be  distributed  according  to  the  Boltzmann  Distribution,  but  the 
particles  spend  the  majority  of  the  time  distributed  around  it.  In  nonequilibrium, 
particles  are  not  distributed  around  the  Boltzmann  Distribution.  The  energy  states 
are  associated  with  the  macroscopic  properties  of  the  system;  the  macroscopic  prop¬ 
erties  are  known  as  long  as  the  temperature  and  the  partition  function  are  known. 
The  partition  function  (Q)  describes  how  the  energy  is  partitioned  across  the  energy 
states,  and  is  defined  as  [12]: 

q  =  E  cie^  <2-5) 

j 

A  particle  can  store  energy  in  different  ways:  as  translational  energy,  rotational  energy, 
vibrational  energy,  or  electronic  energy.  Each  type  of  energy  storage  will  have  its  own 
distributon  function  associated  with  it.  When  particles  go  through  a  shock  wave, 
the  particles’  higher  energy  modes  are  activated,  and  through  collisions  the  particles 
return  to  a  new  equilibrium  state.  The  rate  at  which  particles  return  to  equilibrium 
is  called  the  relaxation  rate  [12],  and  each  type  of  energy  storage  will  have  its  own 
relaxation  rate  dependent  on  types  of  collisions  required  for  reequilibration.  Kinetic 
theory  is  therefore  important  to  discuss  in  order  to  fully  understand  the  behavior  of 
the  gases  through  molecular  collisions. 

There  are  a  few  assumptions  that  go  into  the  derivation  of  the  Boltzmann  equa¬ 
tion  that  must  be  discussed  before  moving  on.  The  density  is  considered  low  enough 


8 


that  only  binary  collisions  occur  and  intermolecular  forces  are  negligible  [13].  In  other 
words,  there  are  so  few  particles  in  the  volume  that  the  likelihood  of  three  molecules 
being  at  the  same  place  at  the  same  time  in  order  to  cause  a  tertiary  collision  is  highly 
unlikely.  Additionally,  since  there  are  so  few  particles,  a  particle’s  path  will  not  be 
altered  simply  by  being  affected  by  another  molecule’s  presence  unless  a  collision 
takes  place.  Another  assumption  is  that  of  molecular  chaos,  which  means  that  the 
velocities  of  two  colliding  particles  are  not  correlated  and  are  independent  of  position. 
Boltzmann  referred  to  the  molecular  chaos  assumption  as  “stosszahlansatz.”  These 
assumptions  are  common  with  DSMC  methodology. 


2. 2  Kinetic  Theory 

Rarefied  conditions  have  been  mentioned  in  the  previous  sections,  but  not  yet 
defined.  The  degree  of  rarefaction  of  a  gas  is  expressed  in  terms  of  the  Knudsen 
number  [4]: 

Kn  =  j  (2.6) 

±j 

where  A  is  the  mean  free  path  is  the  distance  between  molecular  collisions,  and  L  is  a 
characteristic  length  associated  with  the  system.  The  more  rarefied  the  gas  is,  the  less 
likely  it  is  that  particles  will  collide.  The  mean  free  path  will  be  larger,  which  means 
that  particles  will  travel  a  longer  distance  before  colliding,  which  is  the  mechanism 
particles  use  to  equilibrate.  Therefore,  particles  will  not  be  in  equilibrium  with  the 
other  particles  in  the  flowfield  for  longer.  As  a  particle  moves  through  space,  it  sweeps 
out  a  volume  per  unit  time  which  is  ttcPC  [12],  The  number  of  collisions  per  unit 
time  ( 0 )  is  given  by: 

8  =  nd2Cn  (2.7) 

where  d  is  the  diameter  of  the  particle,  C  is  the  average  molecular  velocity.  Using  the 
above  equation,  the  mean  free  path  is: 


A  = 


C 

© 
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Additionally,  the  number  of  collisions  per  volume  (Nc)  is  [4] : 


Nc  =  \n1aTcr 


(2.9) 


Where  ot  is  the  collision  cross-section  and  cr  is  the  relative  velocity  between 
particles.  The  characteristic  length  used  to  find  a  local  Knudsen  number  can  be 
defined  as  [4]: 


L  = 


P 

dp/dx 


(2.10) 


A  high  Knudsen  number  indicates  that  the  mean  free  path  is  near  the  same 
magnitude  as  the  characteristic  length,  which  is  a  condition  that  occurs  at  high  alti¬ 
tudes  when  the  density  is  low  or  if  an  incredibly  small  characteristic  length  is  used, 
such  as  in  nanoscale  simulations. 


Velocity  functions  were  discussed  in  the  preceeding  description  of  the  Boltzmann 
Equation,  but  have  not  yet  been  defined.  The  equilibrium  velocity  distribution  is 
called  the  Maxwellian  distribution: 


m) 


m  \ 
2nkT  ) 


e-^r(ci+c2+ci) 


(2.11) 


where  m  is  the  mass  of  the  particle,  k  is  the  Boltzmann  constant,  T  is  the  tempera¬ 
ture,  and  C  is  the  particle  thermal  velocity  broken  into  its  3  components  [12].  The 
magnitude  of  the  velocity  can  be  investigated  using  the  speed  distribution,  which  is 
derived  from  the  Maxwellian  distribution  [12]: 

^  =  ^{2 Sr)’cV5&Cl  (212) 

The  speed  distribution  (x(C'))  is  written  for  a  gas  with  no  bulk  motion,  which  is  not 
true  for  the  computational  cases  in  this  thesis.  A  method  to  adapt  the  speed  equation 
for  use  in  hypersonic  and  supersonic  flows  has  been  developed,  and  will  be  discussed 
in  Chapter  3.  The  most  probable  thermal  velocity  is  derived  from  the  speed  equation 
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and  is  defined  as: 


The  most  probable  thermal  velocity  (Cmp)  is  the  velocity  that  the  particles  are  most 
likely  to  have  given  the  temperature  of  the  system.  The  Cmp  should  not  be  confused 
with  the  average  velocity,  as  they  are  not  the  same. 

2. 3  Direct  Simulation 

Now  that  the  theoretical  background  that  forms  a  basis  for  DSMC  has  been 
investigated,  DSMC  itself  can  be  discussed.  DSMC  is  a  stochastic  model  that  simu¬ 
lates  a  flowfeld  using  the  Monte  Carlo  statistical  method,  and  is  primarily  used  for 
rarefied  gas  simulations  and  nanoscale  simulations.  The  Monte  Carlo  method  is  used 
to  determine  if  a  collision  occurs.  Once  the  collision  is  accepted,  kinetic  theory  is 
used  to  calculate  the  collision  behavior  of  the  particles.  There  are  a  few  molecular 
models  used  to  describe  particle  collisions  that  have  been  implemented  in  DSMC. 
The  most  basic  is  the  hard  sphere  (HS)  model,  which  is  known  as  the  “billiard  ball” 
model.  In  the  HS  model,  intermolecular  forces  are  neglected,  and  it  is  assumed  that 
molecules  only  interact  with  each  other  by  physically  colliding.  In  other  models,  such 
as  the  Sutherland  model,  that  take  intermolecular  forces  into  account,  molecules  can 
interact  simply  by  being  within  a  certain  radius  of  each  other.  The  Sutherland  model 
adds  weak  atractive  forces  before  the  molecules  contact  [12],  Neither  Sutherlands  nor 
HS  demonstrate  reality.  Molecules  actually  display  an  attractive  force  until  they  get 
within  a  certain  distance  of  each  other,  at  which  point  the  force  becomes  repulsive. 
The  repulsive  force  increases  as  the  distance  between  the  particles  decreases,  and  force 
will  eventually  become  infinitely  strong  [12],  An  improvement  on  the  HS  model  is  the 
variable  hard  sphere  (VHS)  model,  which  is  used  primarily  thoughout  this  project. 
The  VHS  model  varies  the  diameter,  and  therefore  the  collision  cross-section  of  the 
particles  based  on  the  velocity  of  the  particle,  and  will  be  further  discussed  in  another 
section. 
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The  process  of  DSMC  is  fairly  simple,  and  is  based  on  an  assumption  that  the 
motion  and  collision  of  particles  are  decoupled.  The  assumption  allows  the  algorithm 
to  collide  particles  and  then  separately  move  the  particles,  which  will  be  explained  in 
further  detail  when  DSMC  assumptions  are  discussed.  The  DSMC  algorithm  uses  a 
small  time  scale,  usually  on  the  order  of  the  mean  free  time,  which  is  the  time  between 
molecular  collisions.  This  small  time  scale,  usually  on  the  order  of  10-6  seconds,  and 
a  small  cell  size  allow  the  algorithm  to  decouple  the  particle  motion  and  particle 
collision  processes.  The  DSMC  program  progresses  through  a  basic  set  processes  for 
a  specified  number  of  time  steps.  The  process  begins  by  moving  particles  on  the 
velocity  vector  of  each  particle.  The  number  of  collisions  (P)  are  then  calculated 
using  the  probability  equation: 


FNarcrAt 

Vc 


(2.14) 


where  F,y  is  the  ratio  of  actual  particles  per  simulated  particle,  ot  is  the  collision 
cross-section,  cr  is  the  relative  velocity,  At  is  the  change  in  time,  and  Vc  is  the  volume 
of  the  computational  space,  which  for  the  purposes  of  this  equation  is  the  cell  volume. 
Next,  a  randomly  selected  pair  of  particles  from  a  cell  are  chosen  and  the  following 
ratio  is  calculated  and  compared  to  a  randomly  selected  number: 


(Jj^Crp 


(Trcr)r 


>  Ran 


(2.15) 


If  the  randomly  selected  number  is  larger  than  the  ratio,  the  particles  do  not 
collide.  If  the  randomly  selected  number  is  smaller  than  the  ratio,  the  particles 
do  collide.  The  pair  selection  and  comparison  repeats  for  the  calculated  number  of 
collisions  in  the  cell  for  all  cells  in  the  domain.  All  the  particles  are  then  moved  based 
on  their  new  velocity  vectors  (if  collided)  or  the  same  velocity  vector  (if  not  collided) 
and  the  process  is  repeated  until  the  flowheld  reaches  a  steady  state. 
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A  few  notes  to  make  about  the  above  process  before  moving  on:  the  cell  size 
must  be  on  the  order  of  the  mean  free  path  or  smaller  in  order  to  properly  decouple 
the  movement  and  collision  process,  which  can  be  verified  by  making  sure  that  each 
particle  only  collides  once  per  step.  If  a  cell  size  is  on  the  order  of  the  mean  free 
path  and  the  time  step  is  on  the  order  of  the  mean  free  time,  then  a  particle  should 
only  collide  once  per  time  step.  If  the  particles  do  collide  more  than  once  per  step, 
the  assumption  of  decoupling  is  no  longer  valid,  and  either  the  cell  size  or  time  step 
will  need  to  be  reduced.  In  particular,  one  should  look  at  particles  in  the  stagnation 
region  and  shock  layer  (if  one  exists  for  the  particular  flowfield)  since  that  will  be  the 
most  dense  area  of  the  flowfield  and  is  the  limiting  factor  when  designing  a  simulation. 
The  ratio  of  simulated  particles  to  real  particles  can  also  be  increased  to  avoid  more 
than  one  collision  per  particle  per  time  step.  Additionally,  because  of  the  nature  of 
this  statistical  simulation,  if  there  are  not  enough  simulated  particles  in  a  cell  the 
variance  can  be  quite  large.  If  a  flowfield  is  too  rarified,  the  time  steps  or  cell  size  can 
be  increased,  or  the  ratio  of  simulated  particles  to  real  particles  can  be  decreased. 

Bird  made  a  few  key  assumptions  when  he  developed  DSMC  [4],  The  first  is 
local  equilibrium  in  each  cell.  A  flow  through  a  shock  is  not  in  equilibrium,  and  it 
takes  a  number  of  collisions  after  the  shock  to  requilibrate.  Additionally,  flows  within 
a  boundary  layer  are  not  in  equilibrium.  When  there  is  a  gradient  present  in  the 
flowfield,  there  will  be  areas  of  nonequilibrium.  A  shock  layer  is  an  area  between  two 
different  equilibrium  states,  and  there  will  be  a  gradient  present  between  the  states. 
Boundary  layers  exist  because  the  velocity  at  a  wall  is  zero,  so  a  gradient  between  the 
wall  and  the  freestream  develops.  The  flows  in  this  research  contain  both  shock  layers 
and  boundary  layers  and  therefore  are  at  nonequilibrium.  Bird  makes  the  assumption 
that  if  a  cell  size  is  roughly  the  mean  free  path  or  smaller,  the  particles  within  each  cell 
are  in  local  thermodynamic  equilibrium.  Bird  uses  the  following  equation  to  calculate 
the  mean  free  path  using  the  collision  frequency  [4] : 

x-l 

v 
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Where  d  is  the  mean  thermal  speed  of  the  particle,  v  is  the  collision  frequency,  and  n 
is  the  number  density.  Essentially,  there  are  so  few  particles  in  a  small  enough  space 
that  they  can  collide  with  each  other  and  stay  in  equilibrium  even  if  the  flowfield 
itself  is  not  in  equilibrium.  These  cells  create  a  piecewise  equilibrium  where  the 
particles  within  the  cell  are  in  equilibrium,  but  the  particles  are  not  in  equilibrium 
with  the  particles  in  neighboring  cells.  Therefore,  on  a  macroscale  the  flowfield  is  not 
in  equilibrium  even  though  on  a  microscale  the  particles  are  in  equilibrium  with  the 
particles  in  their  immediate  vicinity.  Another  assumption  made  by  Bird  is  within  a 
cell  particle  collisions  are  not  a  function  of  location.  In  other  words,  if  a  cell  is  small 
enough,  one  does  not  need  to  know  where  the  particle  is  within  the  cell  in  order  to 
decide  if  it  can  collide  with  another  particle  within  the  cell.  Since  the  collisions  within 
the  cells  are  not  a  function  of  location,  Boltzmann’s  assumption  of  molecular  chaos 
is  satisfied.  The  assumption  of  molecular  chaos  simply  means  colliding  particles  are 
independent  of  position.  Additionally,  the  two  random  selections  during  the  collision 
process  satisfy  the  molecular  chaos  assumption  that  is  required  with  the  Boltzmann 
equation  [4] .  Another  significant  assumption  of  DSMC  is  that  every  simulated  particle 
represents  many  real  particles,  usually  around  1013,  which  not  only  allows  DSMC  to 
be  computationally  tractable  but  also  provides  the  sampling  pool  of  particles  required 
for  the  statistical  model.  Given  a  small  enough  cell  size,  the  first  two  assumptions 
are  reasonable.  In  conditions  such  as  a  high  Mach  number  flow  over  a  body,  the 
assumption  that  so  many  particles  have  the  exact  same  characteristics  may  not  be 
an  accurate  assumption  [14].  This  last  assumption  will  be  explored  further  as  the 
results  of  the  project  are  discussed  in  Chapter  4. 

2-4  Variable  Hard  Sphere 

As  discussed  previously,  intermolecular  forces  are  neglected  in  the  HS  model, 
and  the  model  is  not  realistic.  The  inverse  power  law  accounts  for  the  change  in 
intermolecular  forces  as  a  function  of  distance  between  particles,  but  it  is  also  com¬ 
putationally  expensive  [15, 16].  The  VHS  model  is  a  combination  of  the  inverse  power 


14 


law  and  HS  models.  The  VHS  model  uses  an  inverse  function  to  relate  the  tempera¬ 
ture  and  the  average  cross-section  of  the  particle.  Therefore,  just  as  the  name  implies, 
a  hard  sphere  is  still  modeled,  but  the  sphere  changes  cross-section  according  to  the 
relationship  [15]: 

aocr(  (2.17) 

where  T  is  temperature  and  (  is  an  exponent  that  is  unique  to  the  particular  gas 
being  examined  and  is  related  to  the  viscosity  index  by: 

lo  =  1/2  +  (  (2.18) 


where  u  is  the  viscosity  index.  The  viscosity  index  for  the  gases  used  in  this  project 
can  be  found  in  Bird  [4],  Equation  2.17  shows  that  as  temperature  increases,  the 
collision  cross-section  will  increase  on  average.  Bird  uses  a  reference  diameter  at  a 
corresponding  reference  temperature  to  calculate  the  diameter  of  the  particle: 


d 


dref 


(2.19) 


where  dref  has  been  defined  in  multiple  references  [4,16],  and  mr  is  the  reduced  mass, 
which  is  defined  as: 

mAmB 

mr  = - 

mA  +  mB 

where  mA  and  mB  are  masses  of  two  colliding  particles.  The  gamma  function  from 
Equation  2.19  is  defined  as: 

OO 

T(j)  =  J  af-'e^dx  (2.21) 

o 

As  one  can  see  in  Equation  2.19,  the  diameter  is  dependent  on  the  relative  velocity, 
and  the  other  values  should  remain  fixed  for  a  given  particle.  As  a  particle’s  velocity 
increases,  the  diameter  (and  therefore  collision  cross-section)  also  increase.  The  SAR 
procedure,  which  will  be  fully  defined  in  the  next  chapter,  also  has  a  velocity  depen- 


(2.20) 
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dency.  Bird  has  proven  that  the  VHS  model  applied  in  a  DSMC  simulation  compares 
well  to  experiment  with  a  normal  shock  wave  with  M=2,  but  at  higher  Mach  numbers 
DSMC  does  not  compare  as  well  to  experimental  data  [11].  Therefore  a  question  that 
must  be  asked  is  whether  SAR  could  be  used  rather  than  or  in  addition  to  VHS  to 
give  accurate  results  for  higher  Mach  flows. 


2. 5  Experimental  Data 

Two  experiments  are  used  for  comparison  in  this  project.  The  first  is  a  1- 
dimensional  normal  shock,  and  the  second  is  the  2-dimensional  axisymmetric  hollow 
cylinder  with  a  step.  These  experiments  are  detailed  in  this  section,  and  in  the 
following  chapter  the  simulations  used  for  comparison  will  be  described. 


2.5.1  1- Dimensional  Normal  Shock.  Alsmeyer  performed  experiments  mea¬ 

suring  the  density  distribution  of  Argon  and  Nitrogen  across  a  normal  shock  wave 
for  Mach  numbers  between  1.55  and  9  [1].  A  stainless-steel  shock  tube  was  used 
for  the  experiment  with  an  inner  diameter  of  150  mm.  The  density  measurements 
were  taken  by  measuring  the  attenuation  of  an  electron  beam  which  was  generated 
by  a  commercial  electron  gun  and  collected  by  a  Faraday  cage.  Once  the  data  was 
collected,  Alsmeyer  plotted  the  normalized  density  profile  across  the  shock  and  the 
inverse  shock  thickness.  The  density  was  normalized  by  [1]: 


Pn 


P~  Pi 
P2  ~  Pi 


(2.22) 


where  p  is  the  measured  density  at  a  point,  pi  is  the  density  before  the  shock,  and  p2 
is  the  density  after  the  shock.  The  normalization  allows  for  comparison  to  Alsmeyer’s 
data  regardless  of  the  initial  density,  which  means  that  many  different  sets  of  exper¬ 
imental  data  can  be  evaluated  with  each  other.  The  error  of  the  normalized  density 
is  accurate  within  ±1%,  and  the  inverse  shock  thickness  is  believed  to  be  accurate  to 
±4%. 
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Figure  2.1:  Inverse  Shock  Thickness  with  Argon  [1] 


Figure  2.1  shows  a  range  of  inverse  shock  thicknesses  at  each  Mach  number, 
which  were  created  by  nine  different  researchers.  Alsmeyer  compiled  the  experimental 
data  and  used  it  as  a  reference  in  Ref  [1].  The  black  line  is  the  curve  fit  of  Alsmeyer’s 
data,  and  the  blue  curve  fit  is  a  6th  order  curve  fit  of  all  the  experimental  data  on 
the  figure.  The  computational  results  will  be  comapred  to  the  blue  curve  fit  of  all  the 
data.  The  data  shows  that  shock  thickness  initially  thins  quickly  as  Mach  number 
increases  until  just  past  Mach  3,  at  which  point  the  shock  thickens  at  a  much  slower 
rate. 


2.5.2  2-Dimensional  Axisymmetric  Hollow  Cylinder.  Davis’  experiment 
that  is  serving  as  a  comparison  is  a  2d  axisymmetric  hollow  cylinder  with  a  step 
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[17].  The  experiments  used  a  heated  nitrogen  wind  tunnel  at  the  Imperial  College 


in  London  in  the  1970s  [18].  The  flow  travels  axially  at  approximately  Mach=21 
with  a  static  temperature  of  23K.  Below  is  a  schematic  of  the  cylinder  used  in  the 
experiment: 


.05m 


Figure  2.2:  Hollow  Cylinder  with  a  Step  [17] 


The  air  flows  axially  along  the  hollow  cylinder  and  impacts  the  6  mm  tall  step 
50  mm  from  the  leading  edge  of  the  cylinder.  The  initial  radius  of  the  cylinder  is  444 
mm,  which  becomes  456  mm  after  the  step.  The  leading  edge  of  the  hollow  cylinder 
is  approximately  0.02  mm  thick  and  has  a  10°  bevel  in  order  to  split  the  flow  around 
the  1  mm  thick  cylinder  [17].  The  cylinder  is  constructed  of  a  copper,  cobalt,  and 
zirconium  alloy  that  has  a  high  thermal  conductivity  in  order  to  maintain  the  uniform 
wall  temperature  of  318K.  The  wall  was  not  water  cooled  for  the  asymmetric  testing 
as  it  was  for  other  studies,  but  the  wall  temperature  was  measured  [17].  The  tem¬ 
perature  measurement  served  as  the  wall  temperature  input  for  the  computational 
study  performed  for  this  thesis.  The  density  values  were  measured  with  an  electron 
beam  probe  [17].  The  nitrogen  gas  fluoresces  when  an  electron  probe  fires  electrons 
through  the  gas,  and  the  intensity  of  the  fluoresence  is  directly  proportional  to  the 
density.  The  flouresence  is  photographed,  and  a  microdensitometer,  which  measures 
the  amount  of  blackness  in  a  photograph,  is  used  to  quantitatively  measure  the  fluo¬ 
resence.  Density  profiles  were  created  at  five  points  along  the  cylinder:  18  mm  before 
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the  step,  3  mm  before  the  step,  at  the  step,  1  mm  after  the  step,  and  7  mm  after  the 
step.  Davis’  density  profiles  can  be  found  in  Figures  2.3  through  2.7. 
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Figure  2.3:  Davis  Density  Profile  at  x=0. 0313m  [17] 


In  Figure  2.3,  the  profile  data  is  taken  ,0187m  upstream  of  the  step.  The  change 
in  density  is  due  to  a  weak  shock  layer  merging  with  the  boundary  layer. 


Figure  2.4:  Davis  Density  Profile  at  x=0. 0462m  [17] 

The  density  profile  in  Figure  2.4  is  taken  0.0038m  upstream  of  the  step.  The 
increased  density  towards  the  wall  is  due  to  the  stagnation  layer.  The  density  peak 
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is  at  approximately  0.032  m  rather  than  0.029  m  in  Figure  2.3,  indicating  a  growing 
boundary  layer  as  distance  from  the  leading  edge  increases. 
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Figure  2.5:  Davis  Density  Profile  at  x=0. 0495m  [17] 


Figure  2.5  shows  the  profile  closest  to  the  stagnation  point  and  therefore  it  has 
the  highest  density  value  at  the  wall.  A  smaller  peak  at  approximately  y=  0.033  m 
shows  the  location  of  the  shock  layer  at  x=0.0495  m. 
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Figure  2.6:  Davis  Density  Profile  at  x=0. 0509m  [17] 
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Past  the  stagnation  point,  the  density  profile  is  similar  to  the  first  profile  at 
x=0.0313  m.  The  change  in  density  is  larger,  as  the  shock  is  stronger  after  the  step. 
However,  both  of  the  profiles  at  x=0.0509  m  and  0.0561  m  show  an  increase  in  density 
close  to  the  wall,  at  approximately  y=0.0305  m. 
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Figure  2.7:  Davis  Density  Profile  at  x=0. 0561m  [17] 


The  shock  is  stronger  6  mm  downstream  from  the  step  compared  to  Figure  2.6, 
as  shown  by  the  change  in  density  at  x=0.0561  m  in  Figure  2.7. 


2. 6  Velocity  Distributions 

The  velocity  distributions  across  a  shock  wave  at  Mach  7.18  were  derived  from 
experimental  data  by  Holtz  and  Muntz  [19].  An  electron  beam  flourescence  technique 
was  used,  just  as  with  the  other  experiments  discussed  in  this  project.  In  order  to 
adequately  resolve  the  regions  in  the  flowfield,  a  Fabry-Perot  etalon  was  implemented 
in  the  observations  [19] .  The  Fabry-Perot  etalon  acts  as  a  filter  and  only  allows  certain 
emission  frequencies  to  be  evaluated.  The  relative  velocity  between  the  sensor  and 
the  molecules  are  calculated  using  the  equation  [19]: 


AA  _  V 
A0  c 


(2.23) 
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Where  AA  is  the  wavelength  shift  due  to  molecular  motion,  A0  is  the  unshifted  wave¬ 
length  of  fluorescence,  V  is  the  relative  velocity,  and  c  is  the  velocity  of  light  [19].  The 
sampling  locations  are  identified  by  a  normalized  density  value  so  that  the  data  may 
be  compared  to  other  experiments  and  computational  results  [19]: 

„  n  —  n\ 

n  = - 

n2  -  ni 

where  n  is  the  normalized  number  density,  n  is  the  number  density  at  the  location, 
ni  is  the  upstream  number  density,  and  n2  is  the  downstream  number  density.  These 
techniques  are  accurate  within  ±4%  [19].  Holtz  and  Muntz  take  data  at  multiple 
points  throughout  the  shock,  but  three  points  were  chosen  to  use  as  a  comparison 
for  this  project:  h  =  .24,  h  =  .333,  and  n  =  .542.  These  points  were  chosen  based 
on  the  changes  occuring  as  the  large  number  of  particles  traveling  at  the  upstream 
velocity  collide  with  other  particles  in  the  shock.  After  h  =  .542,  the  distributions 
take  on  a  fairly  symmetric  appearance  and  do  not  change  significantly.  The  velocity 
distributions  are  normalized  for  ease  of  comparison.  The  velocity  is  divided  by  the 
upstream  bulk  velocity,  and  the  y  axis  is  normalized  by  the  maximum  velocity  function 
so  that  the  distributions  all  peak  at  unity.  The  distributions  chosen  for  comparison  in 
this  project  have  been  reproduced  so  that  they  can  be  plotted  with  the  computational 
results. 


(2.24) 
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Figure  2.8:  Velocity  Distribution  at  n  —  .24  [19] 

The  distribution  in  Figure  2.8  shows  that  most  of  the  particles  are  still  at  the 
bulk  velocity,  which  is  located  at  V jV^.  These  particles  that  are  traveling  at  the  bulk 
velocity  have  not  collided  with  another  particle  in  the  shock  layer,  yet. 


Figure  2.9:  Velocity  Distribution  at  n  =  .333  [19] 

In  Figure  2.9,  the  distribution  is  much  wider  and  there  is  no  longer  a  peak  at 
the  bulk  velocity.  The  majority  of  particles  have  collided  within  the  shock  later  and 
are  traveling  slower  than  the  bulk  velocity.  The  distribution  is  still  slightly  skewed 
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to  the  right  and  not  symmetrical,  indicating  that  there  are  still  a  number  of  particles 
that  have  not  collided  and  are  still  traveling  at  the  bulk  velocity. 


Figure  2.10:  Velocity  Distribution  at  h  =  .542  [19] 

At  the  last  sampling  point  within  the  shock,  the  distribution  is  very  smooth  and 
symmetrical,  indicating  that  at  a  little  over  halfway  through  the  shock  layer  the  vast 
majority  of  particles  have  experienced  a  collision.  The  distributions  before  and  after 
the  shock,  which  are  equilibrium  distributions,  should  also  be  used  for  comparison. 


Figure  2.11:  Velocity  Distribution  Upstream  of  the  Shock  [19] 
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Upstream  of  the  shock,  the  majority  of  the  particles  are  at  the  bulk  velocity, 
and  the  distribution  is  relatively  thin.  Since  the  majority  of  the  particles  are  traveling 
straight  from  left  to  right  at  the  bulk  velocity,  there  are  fewer  collisions,  which  causes 
the  tight  concentration  around  the  bulk  velocity. 


Figure  2.12:  Velocity  Distribution  Downstream  of  the  Shock  [19] 

After  the  shock,  the  vast  majority  of  the  particles  have  collided  and  they  have 
reached  a  new  equilibrium  state.  The  new  bulk  velocity  is  about  half  of  the  upstream 
bulk  velocity.  The  distribution  is  wider,  reflecting  the  fact  that  there  is  a  large  range 
of  particle  velocities  after  the  shock  due  to  collisions. 
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III.  Methodology 

As  stated  in  the  introduction,  the  purpose  of  this  project  is  to  explore  the  applicability 
and  usefulness  of  the  SAR  method.  SAR  will  be  fully  defined  in  this  chapter,  and  the 
computational  experiments  will  be  outlined.  Additionally,  another  variation  of  the 
collision  accept/reject  criteria  will  be  evaluated  called  No  Accept/Reject(NoAR).  As 
the  name  implies,  there  is  no  accept/reject  criteria  in  NoAR:  if  a  pair  of  particles  is 
chosen  to  collide,  then  they  are  automatically  accepted.  A  more  detailed  explanation 
of  NoAR  will  be  found  in  the  next  section. 

The  first  experiment  that  must  be  done  on  these  two  new  algorithms  is  the 
random  walk  test.  Bird  discusses  in  his  book  [4]  that  DSMC  is  a  Markovian  system, 
which  means  that  the  next  time  step  is  only  dependent  on  the  current  time  step. 
A  random  walk  is  a  subset  of  the  Markov  chain,  and  will  not  conserve  molecular 
quantities  such  as  position,  velocity,  and  internal  energies.  Random  walks  can  occur 
in  DSMC  because  the  molecular  quantities  are  conserved  on  the  average  and  because 
of  rounding  that  is  inherent  with  any  computer  program  [4],  Once  SAR  and  NoAR 
have  been  proven  to  contain  no  more  random  walks  than  the  original  DSMC  code, 
the  rest  of  the  experiments  can  begin. 

The  second  experiment  is  the  comparison  to  the  1-dimensional  normal  shock 
experiment  by  Alsmeyer  [1]  using  Bird’s  1-dimensional  code,  called  DSMC1S  [4], 
Alsmeyer  measured  density  and  used  the  measurements  to  calculate  the  inverse  shock 
thickness  of  the  shockwave.  The  inverse  shock  thickness  will  be  calculated  at  several 
Mach  numbers  ranging  from  1.5  to  9,  which  will  not  only  allow  for  comparison  to  ex¬ 
perimental  data,  but  also  provide  insight  into  the  relationship  between  Mach  number 
and  e. 

The  third  is  the  comparison  to  the  hypersonic  flow  over  2-dimensional  axisym- 
metric  cylinder  with  a  step  experiment  that  was  performed  by  Davis  [17].  Davis 
created  density  plots  which  can  be  compared  to  computational  results  using  Bird’s 
DSMC  2-dimensional  axisymmetric  (DSMC2A)  code  [4], 
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The  fourth  comparison  uses  velocity  distributions  created  using  data  from  both 
the  DSMC2A  and  DSMC1S  codes  that  can  be  compared  to  the  velocity  distributions 
from  Holtz  and  Muntz.  Additionally,  velocity  distributions  at  four  different  Mach 
numbers  are  compared  to  each  other  to  understand  the  effects  of  Mach  number  on 
the  velocity  distributions. 

Lastly  HS  simulations  with  SAR  and  NoAR  are  compare  to  the  VHS  simulations 
and  the  experimental  data  for  the  hollow  cylinder  and  the  normal  shock.  Since  both 
SAR  and  VHS  have  Mach  dependencies,  it  was  decided  for  completeness  to  run  HS 
couple  with  SAR. 


3.1  SAR  Method 


The  Smoothed  Accept/Reject  collision  algorithm  varies  only  slightly  from  that 
of  the  original  DSMC  algorithm  by  Bird  [4,  10,  11].  The  collision  accept  criteria 
becomes: 

t~T  mf*  (-  fT  m  f* 

(3.1) 


c JtCt 


±~,  aTCr  >Ran 


(oTcr) 

max  2  (crTCr)r 

where  e  is  a  user  determined  percentage  of  the  ratio  of  collision  cross-section  times 
relative  velocity  over  the  maximum  product  within  the  cell.  Just  as  with  Bird’s 
code,  if  the  random  number  is  larger  than  the  left  side  of  Equation  3.1,  the  collision 
is  rejected,  and  if  the  random  number  is  smaller  the  collision  is  accepted.  Unlike 
Bird’s  code,  however,  random  numbers  that  are  within  the  band  of  the  ratio  plus  or 
minus  |  are  partially  accepted  and  given  a  weighting  between  0  and  1  that  varies 
linearly.  Additionally,  particles  that  are  fully  accepted  are  given  a  weighting  of  1 
and  rejected  particles  are  given  a  weighting  of  zero  [11],  These  weightings  are  used 
to  calculate  the  macroscopic  flowfield  properties.  Only  particles  with  a  non-zero 
weighting  are  considered  when  sampling  the  flowfield,  which  was  hoped  to  make  the 
DSMC  simulations  converge  faster  [11],  using  the  equation: 


Q  = 


E  iWi 


(3.2) 
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where  Qi  is  a  flowfield  parameter  and  wt  is  the  weighting.  Note  that  density  does 
not  use  the  weighting  scheme  in  order  to  conserve  mass.  The  simulations  did  not 
converge  faster,  but  variations  in  the  flowfield  were  noted  for  further  study  [10] .  The 
cases  for  this  project  include  e  values  of  0%,  25%,  50%,  75%,  100%,  and  200%.  The 
e  =  0%  case  has  the  same  accept /reject  criteria  as  the  original  code  and  will  only 
have  weightings  of  zero  and  one  because  there  will  be  no  partially  accepted  collisions. 
The  resulting  flowfield  therefore  should  be  closest  to  the  flowfield  produced  by  Bird’s 
code,  but  not  necessarily  identical  because  of  the  sampling  differences  associated 
with  SAR  compared  to  the  original  DSMC  algorithm.  As  the  e  value  increases,  the 
number  of  rejected  and  fully  accepted  particles  will  decrease  while  the  number  of 
partially  accepted  particles  will  increase  [14] .  In  the  NoAR  code,  if  a  pair  of  particles 
is  selected  to  collide  they  are  accepted  and  given  a  weighting  of  1  [10].  Just  as  with 
SAR,  only  particles  with  a  non-zero  weighting  are  sampled  in  the  NoAR  algorithm, 
and  particles  that  are  not  selected  to  collide  have  a  weighting  of  zero.  The  NoAR 
algorithm  was  developed  because  it  is  the  complete  opposite  of  an  orderly  Markovian 
system  with  an  accept /reject  criteria.  NoAR  is  an  extreme  case  that  higher  e  values 
should  tend  toward.  The  200%  SAR  case  was  selectively  run  to  verify  that  the  results 
do  continue  to  tend  toward  NoAR.  The  expected  results  of  the  study  are  that  Bird 
and  NoAR  produced  flowfield  results  that  are  the  extrema  and  the  SAR  results  are 
in  between,  with  lower  e  values  closer  to  Bird. 

An  illustration  of  the  smoothed  accept/reject  concept  is  seen  below: 

The  blue  line  in  Figure  3.1  represents  the  |  limit  and  the  green  line  is  —  |.  The  pink 
line  is  the  original  Accept/Reject  criteria.  Values  in  Equation  3.1  that  compare  less 
than  the  random  number,  value  are  fully  accepted  and  would  appear  under  the  green 
curve  and  have  a  weighting  of  1.  The  rejected  particles  that  are  above  the  random 
number  would  appear  above  the  blue  line  and  have  a  weighting  of  0.  The  partially 
accepted  particles  would  be  between  the  green  and  blue  lines  and  have  a  weighting 
between  0  and  1. 
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Figure  3.1:  SAR  Illustration 

During  the  previous  work,  a  one-dimensional  shock  simulation  was  compared  to 
experimental  data  using  the  original  algorithm,  SAR,  and  NoAR  [10, 11].  The  results 
of  the  one-dimensional  shock  thickness  experiment  with  Argon  are  shown  in  Figure 
3.2  on  the  next  page. 
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Figure  3.2:  Id  Shock  Thickness  Comparison  with  Experimental  Data 


In  the  Figure  3.2,  inverse  shock  thickness  is  shown  as  a  function  of  Mach  number. 
The  black  line  is  a  curve  fit  of  Alsmeyer’s  experimental  data,  and  the  blue  line  is  the 
curve  fit  for  all  experimental  data  displayed  on  the  graph.  The  figure  also  shows 
inverse  shock  thicknesses  as  calculated  from  the  results  of  the  original  DSMC  code, 
NoAR,  and  e  values  of  10%,  25%,  50%,  75%,  and  90%.  The  calculation  for  inverse 
shock  thickness  is  explained  in  Section  3.3.  The  NoAR  case  results  in  a  thinner  shock 
for  all  cases,  and  the  SAR  results  are  in  between  Bird’s  and  NoAR’s  results.  As  one 
can  see,  the  results  from  Bird’s  original  code  becomes  inaccurate  above  M=3,  when 
Bird’s  code  begins  to  predict  a  thicker  shock  compared  to  experimental  results,  but 
the  SAR  algorithm  is  able  to  match  the  experimental  data.  As  the  Mach  number 
increases,  a  higher  e  value  is  required  to  agree  with  experimental  data,  which  is  a 
result  that  requires  more  study  as  described  in  Section  3.3. 

3.2  Random  Walk  Testing 

In  DSMC,  the  macroscopic  values  are  stored  as  an  average  value  of  all  the  par¬ 
ticles  in  the  cell  [4],  Using  average  values  can  lead  to  random  walks  in  the  Markovian 
system.  A  Markovian  system  is  one  in  which  the  event  at  the  next  step  is  stochasti¬ 
cally  generated  based  on  the  event  at  the  current  time  step  without  information  from 
any  previous  steps  [20].  Random  walks  can  be  tested,  as  Bird  did,  by  simulating  a 
closed  system  where  particles  are  allowed  to  collide  with  each  other  and  the  walls  as 
described  in  10.4  of  Ref  [4],  The  kinetic  energy  of  the  system  is  monitored  for  consis¬ 
tency  over  a  number  of  time  steps.  If  there  is  a  significant  change  in  kinetic  energy, 
random  walks  are  occuring  and  conservation  of  energy  is  not  maintained.  If  kinetic 
energy  remains  nearly  the  same  (one  needs  to  account  for  changes  due  to  round-off 
errors  associated  with  repeating  a  process  many  times),  random  walks  are  not  occur¬ 
ring  in  the  program.  The  particular  test  for  this  project  used  400,000  particles  that 
were  allowed  to  collide  for  1,000  sample  steps.  The  kinetic  energy  and  the  variance  of 
the  kinetic  energy  are  used  for  comparison  between  the  algorithms  and  to  verify  that 
the  random  walks  are  at  a  minimum. 
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3.3  1-Dimensional  Shock  Thickness 


In  addition  to  the  new  work  as  described  later  in  this  chapter,  the  Id  shock 
thickness  comparison  was  reaccomplished,  partially  to  verify  that  the  same  trends  are 
seen  and  partially  to  see  what  additional  information  can  be  garnered  from  the  com¬ 
parison.  The  same  Id  shock  algorithms  used  for  the  speed  distribution  comparisons 
are  used  for  the  shock  thickness  plots.  The  programs  are  used  with  Mach  numbers 
1.5,  2,  3,  6,  and  9.  The  resulting  shock  thickness  for  each  Mach  number  is  then 
graphically  compared  to  Alsmeyer’s  data  [1], 

The  inputs  into  DSMC1S  for  the  project  are  an  initial  temperature  of  293K, 
a  number  density  of  1019,  and  the  ratio  of  simulated  particles  to  real  particles  is 
0.4  *  1016.  The  time  step  is  set  at  0.75  *  10-6  seconds  and  the  cell  size  is  2  mm 
with  a  computational  domain  length  of  0.6  m.  The  project  compares  to  Alsmeyer’s 
argon  experiment,  and  the  properties  of  argon  where  found  in  Bird  [4],  DSMC  uses 
a  velocity  input  with  units  of  meters  per  second  while  Alsmeyer’s  data  is  in  terms  of 
Mach  number.  The  input  velocities  were  calculated  using  the  following  equation: 


M  = 


V 

TWr 


(3.3) 


The  velocities  and  associated  Mach  numbers  are  tabulated  below  for  ease  of  repeati- 
bility: 


Table  3.1:  Velocity  Inputs  Calculated  From  Mach  Number 


Mach 

V  (m/s) 

1.5 

478.25 

2 

637.67 

3 

956.51 

6 

1913.02 

9 

2869.18 

According  to  Vincenti  and  Kruger,  shock  thickness  can  be  related  to  the  pre¬ 
shock  supersonic  velocity  and  the  post-shock  subsonic  velocity  over  the  velocity  gra- 
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dient  [12]: 


(3.4) 


-  'Uq;  lift 

Xn  =  J^{ - 

I  dx  \max 

Where  ua  is  the  before  shock  velocity,  up  is  the  after  shock  velocity,  and  the  denom¬ 
inator  is  the  maximum  velocity  gradient  in  the  shock.  Bird  uses  a  similar  equation, 
but  instead  of  calculating  shock  thickness  in  terms  of  velocity,  he  uses  density  [4], 
Bird’s  equation  is  used  to  find  shock  thickness  in  this  project. 

The  information  needed  for  Bird’s  equation  comes  from  the  data  out  file  that 
is  produced  by  the  DSMC  algorithm,  which  is  then  processed  by  a  MATLAB™code 
written  by  Bentley  and  adapted  by  the  author  to  calculate  the  shock  thickness  [11]. 

3.4  Comparison  to  2-Dimensional  Axisymmetric  Cylinder  Experimen¬ 
tal  Data 

The  third  step  in  the  research  project  is  to  compare  the  results  of  Bird’s  original 
code  to  experimental  data  that  Bird  also  used  for  comparison  [4, 14]  by  Davis.  Sim¬ 
ulations  for  this  portion  of  the  project  were  conducted  using  Bird’s  original  code  [4], 
a  NoAR  case,  and  a  range  of  e  values  which  are  0%,  25%,  50%,  75%,  100%,  and 
200%.  Bird’s  2d  axisymmetric  code,  and  all  the  DSMC  codes  used  in  this  project, 
can  only  create  surfaces  along  a  cell  boundary  and  cells  can  only  be  rectangular  [4], 
Therefore  the  leading  edge  will  not  have  a  bevel  in  the  simulations  which  may  have  a 
minor  impact  on  the  simulation  -  that  impact  will  be  taken  into  consideration  when 
reviewing  the  results.  The  dimensions  of  the  computational  space  are  0.15  m  in  the 
x-directon  and  0.085m  in  the  y-direction  in  order  to  match  the  experiment  and  sub¬ 
sequent  simulations  performed  by  Davis  [17].  The  computational  domain  in  relation 
to  the  experimental  set-up  for  the  hollow  cylinder  is  shown  below  in  Figure  ?? 

The  hollow  cylinder  is  1  mm  from  the  bottom  of  the  computational  domain,  and 
the  leading  edge  begins  at  the  left  side  of  the  domain.  The  cell  size  for  the  simulations 
is  1  mm  in  the  x  and  y  directions.  A  grid  with  cell  dimensions  of  0.5  mm,  which  is 
on  the  order  of  the  mean  free  path,  was  also  used  in  the  study  to  see  if  the  solution 
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85mm 


Figure  3.3:  Computational  Domain 

is  improved.  The  hollow  radius  of  the  cylinder  is  not  part  of  the  computational  space 
in  the  simulation,  but  the  radius  is  not  neglected.  Bird’s  2-dimensional  axisymmetric 
code  accounts  for  an  internal  radius  that  is  not  in  the  computational  space  [4], 

In  addition  to  comparing  the  results  to  experimental  data,  the  results  have 
been  compared  to  each  other.  Contour  plots  generally  provide  a  qualitative  method 
of  comparing  results.  In  order  to  improve  upon  the  comparison,  the  contour  plots  have 
been  generated  by  subtracting  the  result  of  each  of  the  SAR  and  NoAR  algorithms 
from  Bird’s  results,  and  dividing  by  Bird’s  result  and  multiplying  by  100  in  order 
to  plot  the  percentage  of  difference  between  the  results.  For  example,  the  density 
percent  difference  between  Bird  and  the  results  of  SAR  with  e  =  0%  is  found  by: 


PBird  PEps 00 nn  /Q 

P%diff  = - 100  (3.5) 

PBird 

There  are  some  cases  where  the  value  of  density  will  be  zero,  in  particular  where 
the  step  is.  In  Bird’s  code,  the  wall  boundaries  are  defined  in  the  data  subroutine, 
but  these  areas  with  no  flow  are  a  part  of  the  computational  domain,  and  data  is 
written  out  for  these  locations.  Therefore,  the  areas  of  no  flow  need  to  be  taken  into 
account.  For  these  cases,  the  density  is  assumed  to  be  10-6,  which  gives  a  large  value 
for  the  percent  difference.  The  difference  calculations  were  performed  for  density, 
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u-velocity,  and  temperature.  These  contour  plots  also  give  an  idea  of  which  cell  the 
velocity  data  should  be  exported  from  for  the  speed  distributions  in  order  to  see  the 
largest  difference  between  results.  It  would  make  sense  that  the  cells  with  only  small 
differences  in  the  u-velocity  would  have  similar  velocity  distributions  while  cells  with 
larger  differences  will  have  disimilar  velocity  distributions.  Choosing  of  the  sample 
cells  is  further  discussed  in  Section  3.5. 

3.5  Velocity  Distribution  Comparisons 

The  velocity  distributions  as  described  in  Section  2.6  will  provide  further  insight 
into  the  behavior  of  the  particles  as  they  travel  through  the  shock  layer.  Before  the 
shock,  the  majority  of  the  particles  are  traveling  at  the  bulk  velocity,  which  creates 
a  thin  profile.  Within  the  shock,  the  particles  transition  from  traveling  at  the  bulk 
velocity  to  traveling  at  a  lower  speed  via  molecular  collisions.  The  distributions  in 
Section  2.6  show  the  behavior  of  the  particles  as  they  collide  and  eventually  reequili¬ 
brate  to  a  new  flow  velocity.  Velocity  distributions  resulting  from  the  DSMC,  SAR, 
and  NoAR  simulations  will  be  compared  to  Holtz  and  Muntz’s  data  at  Mach  7.18, 
and  velocity  distributions  will  also  be  created  for  Mach  1.5,  3,  6,  and  9  in  order  to 
observe  the  change  in  the  distributions  with  Mach  number.  As  the  Mach  number  in¬ 
creases,  the  shock  becomes  stronger  and  the  perturbation  from  the  equilibrium  state 
becomes  larger.  The  velocity  distributions  should  reflect  the  change  in  the  amount  of 
nonequilibrium,  and  DSMC,  SAR,  and  NoAR  may  handle  the  changes  differently. 

The  distributions  will  all  be  normalized  just  as  Holtz  and  Muntz  normalized 
their  distributions  [19].  Using  the  normalized  distributions  gives  a  height  of  1  for 
all  the  distributions  and  compares  the  particle  velocities  to  the  bulk  velocity  of  the 
flowfield,  which  allows  for  direct  comparisons  regardless  of  the  input  bulk  velocity 
used  in  the  DSMC,  SAR,  or  NoAR  codes.  The  histogram  program  for  the  unweighted 
velocities  uses  the  histogram  function  in  MATLAB™to  calculate  the  distributions, 
which  are  then  normalized.  The  weighted  velocities  required  a  different  approach. 
With  the  unweighted  velocity  histogram,  if  a  particle’s  velocity  falls  within  a  certain 
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bin,  then  one  is  added  to  the  number  of  particles  within  that  bin.  The  weighted 
particles  however  do  not  all  have  a  weighting  of  one.  Therefore,  if  a  weighted  particle’s 
velocity  falls  within  a  bin,  the  weighting  is  added  to  the  number  of  particles  in  the 
bin. 

Speed  distribution  studies  were  initially  accomplished  with  DSMC1S,  and  the 
modifications  to  the  program  with  SAR  and  NoAR.  The  particle  velocity  components 
and  weights  at  locations  within  the  flowfield  were  exported  into  a  file.  The  compu¬ 
tational  domain  has  300  cells,  with  the  shock  occuring  at  the  middle.  The  shock  is 
at  x=0  m,  the  left  boundary  is  at  x=-0.3  m,  and  the  right  boundary  is  at  x=0.3  m. 
Mach  1.5,  3,  6,  and  9  simulations  have  three  sampling  locations:  before  the  shock,  in 
the  shock,  and  after  the  shock.  The  preshock  sampling  location  is  10  cells  into  the 
domain  from  the  left  side  of  the  domain,  at  x=-0.28  m.  The  location  within  the  shock 
was  chosen  at  cell  150  at  x=0m,  and  the  post  shock  location  is  cell  290  at  x=0.28  m. 
The  Mach  7.183  distributions  were  taken  at  sample  locations  based  on  the  location  of 
the  normalized  density.  The  DSMC,  SAR,  and  NoAR  simulations  were  all  run  10,000 
time  steps,  and  the  normalized  density  was  calculated  for  each  case.  The  simulations 
were  then  run  again  and  the  particle  velocities  were  written  at  the  locations  of  the 
normalized  densities  previously  mentioned  in  Section  2.6. 

The  DSMC2A  code  is  also  used  to  investigate  the  velocity  distributions.  The 
flowfield  is  very  complicated,  and  the  velocity  distributions  are  difficult  to  predict. 
The  same  flowfield  complications  make  it  difficult  to  choose  which  cells  should  be 
sampled  for  the  distributions.  Four  locations  were  initially  chosen  as  sample  cells 
for  the  velocity  distributions,  which  are  495,  795,  1255,  and  1295  as  seen  in  Figures 
3.4  and  3.5.  Cell  1295  was  in  the  boundary  layer  for  some  of  the  simulations,  which 
affected  the  comparison  of  velocity  distributions,  so  1255  and  1295  were  replaced 
with  1705  and  1755.  Cells  350  and  3669  were  chosen  by  using  the  percent  difference 
contour  plots  to  determine  the  areas  of  the  flow  with  the  largest  differences.  Not 
entirely  surprisingly,  the  biggest  differences  were  found  at  the  stagnation  point  (350) 
and  in  the  shock  layer  (3669).  The  points  are  identified  by  the  cell  number  that  they 
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are  assigned  in  DSMC.  DSMC  assigns  each  cell  a  number  starting  with  one  in  the 
bottom  left  hand  corner,  and  goes  left  to  right  and  then  up  to  the  next  row  of  cells. 
The  figure  below  shows  the  location  of  the  points  in  reference  to  the  flowfield: 


Figure  3.4:  Sampled  Points  on  Velocity  Contour  Plot 


The  location  of  the  shock  is  more  easily  visible  in  Figure  3.5,  which  is  a  contour 
plot  of  density  rather  than  velocity. 
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Figure  3.5:  Sampled  Points  on  Density  Contour 
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Figures  3.4  and  3.5  show  that  the  recirculation  region  is  represented  by  cells  495 
and  795.  Cell  350  is  in  the  stagnation  region,  and  1255  and  1705  are  just  after  the 
step.  Cells  1295  and  1755  will  show  the  behavior  well  past  the  step,  and  cell  3669  will 
allow  the  investigation  of  the  shock  layer  velocity  distributions.  These  points  were 
sampled  using  Bird’s  code,  NoAR,  and  SAR  with  e  values  of  00%,  25%,  50%,  75%, 
and  100%.  The  distributions  are  compared  visually  and  with  the  statistical  equations 
discussed  earlier  in  the  section. 

After  creating  all  the  plots  and  visually  comparing  them,  it  was  determined  that 
the  histograms  are  only  minutely  different.  A  more  precise  way  to  compare  them  must 
be  developed.  Evaluating  the  characteristics  of  the  distributions,  such  as  maximum 
value,  the  average  velocity,  Cmp,  and  moments  of  the  histogram,  specifically  kurtosis 
and  skewness.  The  maximum  value  of  n  will  be  located  in  the  velocity  bin  with  the 
most  probable  velocity.  If  Cmp  matches  the  location  of  the  maximum  n  value,  then 
the  histogram  programs  are  working  correctly.  A  distribution  with  high  kurtosis  will 
be  tall,  thin,  and  pointed  at  the  peak,  while  a  low  kurtosis  distribution  will  be  short, 
and  broad.  Sample  kurtosis  is  calculated  by  [20]: 


-4_  ITJU  (au-T)4 

‘  (iE?=i(*i-s)2)5 


Skewness  is  the  measure  of  symmetry  around  the  mean.  A  distribution  with  no 
skewness  will  be  perfectly  symmetric  about  the  mean.  A  distribution  with  negative 
skewness  is  asymmetric,  favoring  the  left  of  the  mean  and  a  distribution  with  a  positive 
skewness  favors  the  right  side  of  the  mean.  Skewness  of  a  sample  distribution  can  be 


calculated  by  [20]: 


(3.7) 


3.6  Hard  Sphere  Comparison 

In  order  to  verify  if  the  SAR  method  can  be  used  as  a  Mach  dependent  param¬ 
eter,  rather  than  VHS,  that  can  be  varied  to  produce  a  more  accurate  result,  it  is 
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necessary  to  use  the  HS  model  with  SAR  and  compare  to  Bird’s  original  code  for  a 
given  Mach  number.  The  DSMC1S  and  DSMC2A  codes  were  both  altered  for  HS. 
The  collision  cross-section  of  the  particle  is  calculated  using  [4] : 

ot  =  7rd2  (3.8) 

The  diameter  is  constant,  unlike  the  diameter  in  Equation  2.19.  The  results  of  the 
hard  sphere  model  cases  will  be  compared  both  to  Bird’s  original  code  with  VHS 
and  SAR  and  NoAR  codes  with  VHS.  For  DSMC2A,  HS  will  be  compared  using  the 
density  plots  from  Davis.  The  DSMC1S,  HS  codes  will  be  evaluated  using  line  plots 
and  inverse  shock  thickness  calculations. 
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IV.  Results 


4-1  Random  Walk  Testing 

The  random  walk  testing  results  include  variance  and  average  kinetic  energy.  If 
the  average  kinetic  energy  stays  relatively  constant,  there  are  no  random  walks  in  the 
system  and  the  algorithms  are  conserving  energy. 
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Figure  4.1:  Variance  for  Random  Walk  Test 

As  can  be  seen  in  Figure  4.1,  the  variance  for  all  three  cases  are  on  the  same 
order  of  magnitude.  The  average  kinetic  energy  stays  constant  at  about  0.0004  while 
the  variance  is  between  —  6*  10-12  and  6*  10-12.  The  small  amount  of  variance  is  most 
likely  caused  by  rounding  errors  due  to  the  computations  performed  by  the  computer, 
rather  than  problems  with  the  codes. 
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4-2  1- Dimensional  Shock  Experimental  Comparisons 

4-2.1  Inverse  Shock  Thickness.  A  few  small  changes  were  made  in  the 
sampling  subroutine  of  DSMC1S  from  previous  research,  so  the  first  part  of  the  inves¬ 
tigation  was  to  verify  that  the  inverse  shock  thickness  plot  is  still  showing  the  same 
behavior  as  was  noted  before. 

Just  as  with  Figure  3.2,  the  black  line  in  Figure  4.2  is  Alsmeyer’s  curve  fit,  and 
the  blue  line  is  the  6th  order  curve  fit  for  all  the  experimental  data.  The  crosses  all 
represent  simulation  results,  with  the  red  being  DSMC,  blue  is  SAR  00%,  pink  is  SAR 
25%,  yellow  is  SAR  50%,  cyan  is  SAR  75%,  green  is  SAR  100%,  and  black  is  NoAR. 
Through  Mach  3,  Bird’s  code  matches  well  with  the  blue  curve  fit,  but  at  Mach  4, 
Bird  begins  to  over  predict  the  shock  thickness.  Bird’s  symbol  cannot  always  be  seen 
on  the  figure  because  the  SAR  00%  data  is  colocated. 

An  effort  was  then  made  to  curve  fit  the  e  input  values  that  would  provide  the 
correct  inverse  shock  thickness  for  each  given  Mach  number,  and  the  resulting  plot 
can  be  seen  in  Figure  4.3.  The  initial  value  of  e  at  each  Mach  number  was  found  by 
interpolating  the  e  values  using  Figure  4.2.  These  values  were  then  used  as  an  input  in 
a  SAR  simulation,  and  the  output  was  used  to  calculate  the  inverse  shock  thickness. 
The  resulting  inverse  shock  thickness  was  then  evaluated,  and  a  new  e  value  was 
chosen  if  the  inverse  shock  thickness  did  not  match  the  blue  curve  fit.  The  process 
was  repeated  until  the  inverse  shock  thickness  for  each  Mach  number  all  matched  the 
blue  curve  fit  line  in  Figure  4.2.  Figure  4.3  is  the  result  of  the  e  fitting  process  just 
described. 
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Figure  4.2:  Inverse  Shock  Thickness  vs  Mach  Number  [1] 
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Figure  4.3:  Inverse  Shock  Thickness  vs  Mach  Number  SAR  Fit 


The  cyan  crosses  in  Figure  4.3  are  the  calculated  results.  The  Mach  number 
and  their  associated  e  values  are  displayed  in  the  Table  4.1.  The  additional  Mach 
numbers  were  added  in  areas  where  the  change  in  e  versus  Mach  number  were  large 
in  order  to  get  a  good  understanding  of  the  behavior  in  high  gradient  portions  of  the 
curve.  Table  4.1  shows  the  values  plotted  in  Figure  4.3. 

Table  4.1:  e  Values  Fit  to  Mach  Number 


Mach 

*(%) 

1.5 

0 

2 

0 

2.5 

.5 

3 

5 

3.5 

10 

4 

15 

4.5 

18 

5 

22 

6 

27 

7 

28 

8 

30 

8.25 

30.5 

8.5 

33 

8.75 

34 

9 

45 

Since  Bird’s  results  fit  through  Mach  3,  the  first  two  points  in  the  curve  fit  are 
zero,  which  will  negatively  affect  the  curve  fit.  Two  curve  fits  were  therefore  created, 
one  that  includes  the  first  5  points,  and  one  that  does  not.  The  two  curves  in  Figure 
4.4  allow  one  to  tell  which  curve  fit  will  work  best. 

The  orange  line  in  Figure  4.4  represents  the  curve  fit  over  the  entire  data  set, 
the  green  line  represents  the  curve  fit  of  data  above  Mach  4,  and  the  blue  crosses  are 
the  points  used  to  calculate  the  curve  fits.  One  can  see  that  the  curve  fit  for  the  whole 
data  set  does  not  match  the  points  between  Mach  8  and  Mach  9,  but  the  curve  fit 
for  just  the  higher  Mach  data  does  match.  Note  the  significant  increase  in  e  required 
to  maintain  the  correct  inverse  shock  thickness  at  higher  Mach  numbers.  In  order 
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Mach  Number 


Figure  4.4:  Inverse  Shock  Thickness  vs  Mach  Number  Curve  Fit 

to  better  understand  the  relationship  between  the  algorithms  and  Mach  number  in 
producing  the  correct  inverse  shock  thickness,  a  graph  was  created  that  shows  the 
percent  difference  between  the  experimental  inverse  shock  thickness  and  the  inverse 
shock  thickness  calculated  by  the  algorithms.  The  experimental  values  used  for  the 
graph  came  from  the  blue  curve  fit  of  all  the  experimental  data  on  Figure  4.2. 

In  Figure  4.5,  Bird’s  and  e  =  0%’s  results  match  the  experimental  data  at 
Mach  1.5,  but  then  underpredict  the  shock  thickness  at  Mach  2.  After  Mach  2,  the 
two  algorithms  consistently  predict  a  thicker  shock  compared  to  experimental  data, 
e  =  25%  underpredicts  the  inverse  shock  thickness  through  Mach  5,  but  at  Mach 
6,  is  very  close  to  the  experimental  data,  which  matches  Figure  4.4.  The  higher  e 
values  and  NoAR  predict  a  thinner  shock  at  every  Mach  number  in  Figure  4.5.  The 
overall  behavior  of  Figure  4.5  matches  Figure  4.4.  The  percent  differences  steadily 
rise  until  Mach  6  in  Figure  4.5,  and  flattens  out  until  Mach  8,  at  which  point  the 
percent  difference  increases  drastically.  The  curve  fit  in  Figure  4.4  increases  also  until 
Mach  6,  becomes  relatively  flat  until  Mach  8,  then  spikes  at  Mach  9.  One  would 
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Figure  4.5:  Inverse  Shock  Thickness  Percent  Difference  vs  Mach  Number 

expect  that  these  two  graphs  would  match,  since  both  of  these  graphs  are  different 
ways  to  show  the  relationship  between  Bird,  SAR  with  different  e  input  values,  and 
NoAR  to  an  inverse  shock  thickness  associated  with  Mach  numbers.  The  relationship 
should  be  extended  to  higher  Mach  number  data  set  to  best  confirm  the  relationship 
or  modify  it  as  needed  for  higher  Mach  numbers. 

Now  that  it  has  been  established  that  there  is  a  Mach  dependency,  it  is  necessary 
to  find  out  what  causes  the  dependency.  The  first  thing  to  look  at  is  the  collision 
ratio,  which  is  the  ratio  of  pairs  of  particles  accepted  (fully  and  partially)  over  the 
total  number  of  pairs  selected  as  possible  collision  partners. 

In  Figure  4.6,  the  blue  diamond  is  e  =  100%,  the  red  square  is  e  =  75%,  the 
green  triangle  is  e  =  50%,  the  purple  x  is  e  =  25%,  the  blue  start  is  e  =  0%,  and 
the  orange  circle  is  Bird.  NoAR  by  definition  always  has  a  collisions  ratio  of  1.  Bird 
has  the  smallest  collision  ratio,  which  is  expected  since  it  has  a  binary  accept/reject 
criteria,  e  =  0%  provides  the  same  values  for  the  collision  ratio  as  Bird.  As  e  increases, 
so  does  the  collision  ratio.  The  e  value  will  increase  the  number  of  partially  accepted 
collisions,  which  increases  the  overall  collision  ratio.  As  Mach  number  increases,  the 
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Figure  4.6:  Collision  Ratio  vs  Mach  Number 

collision  ratio  decreases,  and  begins  to  reach  an  asymptotic  limit.  This  limit  is  due 
to  the  change  of  the  velocity  distribution  with  velocity. 


Figure  4.7:  Change  in  Velocity  Distribution  with  Mach  Number 

In  Figure  4.7,  the  black  line  is  Mach  1.5,  the  magenta  line  is  Mach  3,  the  red 
line  is  Mach  6,  the  blue  line  is  Mach  9,  and  the  cyan  line  is  Mach  12.  The  dotted  black 
and  cyan  lines  represent  the  upper  limit  of  the  SAR  criteria,  and  the  dashed  line  is 
the  lower  limit  of  the  SAR  criteria  all  for  an  e  =  25%  input  value.  One  can  see  how 
at  lower  velocities  the  distribution  is  tall  and  thin,  but  as  the  velocity  increases  the 
distribution  becomes  shorter  and  wider.  The  change  between  Mach  1.5  and  Mach  3 
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is  drastic,  but  as  the  velocities  become  larger,  the  change  in  the  distribution  becomes 
smaller.  These  smaller  changes  result  in  smaller  differences  in  collision  ratio  as  the 
Mach  number  increases.  Additionally,  as  the  velocity  increases,  the  maximum  gtct 
value  is  going  to  increase,  which  will  decrease  the  likelihood  of  a  particle  pair  having 
a  large  ratio  in  Equation  2.15,  which  will  decrease  the  number  of  collisions  that  are 
accepted.  This  behavior  is  small  though,  and  does  not  explain  the  drastic  change 
in  e  as  Mach  number  increases  required  to  match  experimental  data.  The  velocity 
distributions  will  allow  for  further  investigation  of  the  behavior  of  DSMC,  SAR,  and 
NoAR. 

4-2.2  Velocity  Distribution  Comparisons.  The  velocity  distribution  compar¬ 
isons  provide  insight  into  the  differences  of  the  algorithms,  especially  with  how  well 
the  algorithms  keep  with  the  assumption  of  local  equilibrium.  As  previously  stated, 
local  equilibrium  within  the  computational  cells  can  be  assumed  if  the  cell  size  in 
on  the  order  or  smaller  than  the  mean  free  path.  The  cell  size  for  the  1-dimensional 
algorithms  is  one  order  of  magnitude  smaller  than  the  mean  free  path.  The  Mach  7.18 
case  can  be  compared  to  Holtz  and  Muntz’s  experimental  data.  These  distributions 
are  taken  using  the  velocity  traveling  in  the  x-axis  in  order  to  compare  appropriately 
to  the  experimental  data. 

In  Figure  4.8,  Bird’s  result  is  the  tan  color  with  the  x  symbol,  e  =  0%  is  gray 
with  the  line,  the  e  =  50%  is  the  red  line  with  the  triangle,  e  =  100%  is  the  green 
dashed  line,  NoAR  is  the  blue  line  with  the  upside  down  triangle,  and  Holtz  and 
Muntz’s  data  is  the  blue  circle.  All  of  the  distributions  in  Figure  4.8  are  unweighted. 
The  SAR  distributions  are  created  by  using  the  modified  accept/reject  criteria,  but 
the  weightings  are  neglected,  which  allows  for  an  understanding  of  the  effect  of  the 
accept/reject  portion  of  the  algorithm  by  itself. 
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Unweighted  Upstream  Velocity  in  x-direction 


Figure  4.8:  Unweighted  Velocity  Distributions  Upstream  of  Shock 


The  distributions  in  Figure  4.8  are  taken  upstream  of  the  shock,  and  at  this  point 
the  distributions  are  thin.  These  cells  are  in  equilibrium,  and  all  of  the  distributions 
follow  Holtz  and  Muntz’s  experimental  data  similarly  well. 


Weighted  Upstream  Velocity 


Figure  4.9:  Weighted  Velocity  Distributions  Upstream  of  Shock 
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Figure  4.9,  the  distributions  are  created  using  the  particles  weights  based  on 
the  SAR  methodology.  Just  as  with  the  previous  figure,  all  of  the  distributions  match 
the  experimental  data  well. 


Unweighted  Velocity  at  n=.24  in  x-direction 


Figure  4.10:  Unweighted  Velocity  Distributions  at  h  =  0.24 

Figure  4.10  shows  the  unweighted  distributions  about  a  quarter  of  the  way 
through  the  shock  structure.  Bird  and  the  unweighted  SAR  distributions  severely 
underestimate  the  number  of  collided  particles.  There  is  a  large  number  of  particles 
at  the  bulk  velocity,  and  few  particles  traveling  slower  than  the  bulk  velocity.  The 
e  =  0%  shows  more  particles  traveling  slower  than  the  bulk  velocity,  and  as  the  e  value 
increases,  the  number  of  particles  tends  towards  NoAR.  The  NoAR  distribution  has 
the  lowest  number  of  particles  traveling  slower  than  the  bulk  velocity,  indicating  that 
the  number  of  collided  particles  is  fewer.  The  changes  in  accept/reject  criteria  allow 
for  a  wider  range  of  particles  to  collide  that  would  have  otherwise  been  unable  to, 
which  causes  minor  changes  in  the  distribution,  as  seen  in  the  area  of  the  distribution 
to  the  left  of  the  peak  at  unity. 

Figure  4.11  shows  a  similar  peak  as  the  unweighted  distributions,  but  the  por¬ 
tion  of  the  distributions  to  the  left  of  the  peak  are  very  different  from  the  unweighted 
cases.  The  weighted  e  =  0%  distribution  matches  well  to  the  experimental,  but  pre- 
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Weighted  Velocity  at  n=.24 


Figure  4.11:  Weighted  Velocity  Distributions  at  h  =  0.24 

diets  slightly  more  particles  at  the  lower  velocities  than  the  experimental  observations. 
The  e  =  50%  distribution  underpredicts  the  number  of  particles  traveling  less  than 
the  bulk  velocity,  as  do  e  =  100%  and  NoAR.  As  e  increases,  the  number  of  parti¬ 
cles  traveling  less  than  the  bulk  velocity  decreases,  which  is  most  likely  due  to  the 
weighting.  As  the  band  of  partially  accepted  collisions  increases,  more  particles  are 
given  a  weighting  between  zero  and  one,  and  these  partial  weightings  will  decrease 
the  number  of  particles  compared  to  e  =  0%. 

In  Figure  4.12,  the  velocity  distributions  are  taken  at  approximately  one  third 
of  the  way  through  the  shock  layer.  The  experimental  plot  widens  as  more  particles 
have  now  collided  and  are  traveling  at  velocities  other  than  the  bulk  velocity.  The 
unweighted  distributions  still  show  a  thin  peak  at  the  bulk  velocity  and  underpredicts 
the  number  of  particles  that  are  going  slower  than  the  bulk  velocity.  The  number  of 
particles  to  the  left  of  the  peak  has  increased  compared  to  Figure  4.10,  indicating  that 
more  particles  have  collided  as  the  particles  have  move  farther  into  the  shock  layer. 
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Unweighted  Velocity  at  n=.333  in  x-direction 


Figure  4.12:  Unweighted  Velocity  Distributions  at  h  =  0.333 


Weighted  Velocity  at  n=.333 


Figure  4.13:  Weighted  Velocity  Distributions  at  h  =  0.333 


The  weighted  distributions  in  Figure  4.13  again  show  a  much  better  comparison 
to  experimental  data  than  the  unweighted  distributions.  All  of  the  distributions 
underpredict  the  number  of  particles  to  the  left  of  the  peak,  but  match  the  peak  well. 
As  seen  previously,  the  e  =  0%  case  underpredicts  the  least,  and  NoAR  underpredicts 
the  most. 
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Unweighted  Velocity  at  n=.542  in  x-direction 


Figure  4.14:  Unweighted  Velocity  Distributions  at  h  =  0.542 


The  unweighted  distributions  in  Figure  4.14  match  the  experimental  data  very 
well  to  the  left  of  the  peak  at  one,  and  the  width  of  the  distributions  compare  well  to 
the  experimental  data.  Bird  and  e  =  0%  have  the  most  pronounced  peak  at  one,  and 
as  e  increases,  the  peak  is  reduced.  NoAR  does  not  show  a  peak,  but  the  slope  does 
not  quite  match  the  experimental  data. 


Weighted  Velocity  at  n=.542 


Figure  4.15:  Weighted  Velocity  Distributions  at  n  —  0.542 
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Figure  4.15  shows  that  all  of  the  weighted  distributions  match  very  well  to  the 
experimental  data.  The  experimental  velocity  distribution  is  very  near  to  the  post 
shock  equilibrium  distribution,  and  as  will  be  seen  in  the  next  figure,  there  is  little 
change  after  this  point  in  the  shock  structure.  All  of  the  weighted  distributions  show 
that  enough  collisions  have  occured  in  order  to  compare  favorably  to  the  experimental 
data.  Even  the  unweighted  distributions  shown  in  Figure  4.14  are  fairly  close  to 
matching  the  experimental  distribution. 


Unweighted  Downstream  Velocity  in  x-direction 


Figure  4.16:  Unweighted  Velocity  Distributions  Downstream  of  Shock 

Downstream  of  the  shock,  in  Figure  4.16,  the  flowheld  has  returned  to  an  equi¬ 
librium  state.  The  unweighted  distributions  again  are  very  similar  to  the  experimental 
data,  which  confirms  the  idea  that  in  equilibrium  the  algorithms  produce  similar  re¬ 
sults  because  the  collision  integral  is  not  a  factor. 

In  Figure  4.17,  there  is  very  little  difference  between  these  distributions,  and 
they  all  compare  well  to  the  experimental  data.  Just  as  with  the  unweighted  distri¬ 
butions  downstream  of  the  shock,  now  that  the  flow  has  returned  to  an  equilibrium 
state,  all  of  the  algorithms  produce  comparable  results. 
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Weighted  Downstream  Velocity 


Figure  4.17:  Weighted  Velocity  Distributions  Downstream  of  Shock 

The  results  of  the  experimental  comparison  show  that  in  equilibrium,  the  results 
of  the  algorithms  are  very  similar,  but  in  nonequilibrium  within  the  shock  layer, 
clear  differences  can  be  seen.  Bird  and  the  unweighted  SAR  and  NoAR  distributions 
show  that  the  algorithms  do  not  predict  the  collision  rate  accurately,  which  results 
in  particles  maintaining  the  bulk  velocity  longer  through  the  shock  than  is  seen  in 
experiment.  The  weighted  SAR  distributions  better  match  the  experimental  data, 
with  the  lower  e  =  0%  matching  best  compared  to  e  =  100%  and  NoAR.  NoAR 
matches  the  worst  of  the  weighted  distributions  through  the  shock,  even  though  NoAR 
has  the  highest  collision  rate. 

4-3  1 -Dimensional  Data 

4-3.1  Line  Plots.  Line  plots  were  created  using  Mach  1.2  and  Mach  9  cases 
for  density,  u-velocity,  and  temperature  in  order  to  gain  an  understanding  of  how 
the  SAR  algorithm  affects  the  flowfield  parameters  at  the  lowest  and  highest  Mach 
numbers  evaluated  in  this  study.  Mach  1.2  was  the  first  case  evaluated: 

In  Figure  4.18,  Bird’s  result  is  green,  e  =  0%  is  red,  e  =  25%  is  blue,  e  =  50% 
is  brown,  e  =  75%  is  light  blue,  e  =  100%  is  orange,  and  NoAR  is  purple.  The  color 
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Figure  4.18:  Density  Line  Plot  for  Mach  1.2  Normal  Shock 

scheme  for  each  case  is  kept  throughout  this  report  for  ease  of  reading.  Bird’s  results 
and  e  =  00%  do  not  he  on  top  of  each  other.  However,  Bird  and  NoAR  are  still  the 
extrema  with  the  SAR  results  in  between.  The  same  is  true  for  the  velocity  in  Figure 
4.19. 


Figure  4.19:  U-Velocity  Line  Plot  for  Mach  1.2  Normal  Shock 
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Magnifying  the  figure  at  the  top  of  the  profile,  it  can  be  seen  that  all  the  results 
of  e  =  75%  and  higher  lie  on  the  NoAR  results.  Additionally,  the  NoAR  result  does 
not  show  as  much  variance  as  the  other  cases. 


Figure  4.20:  U-Velocity  Line  Plot  for  Mach  f.2  Normal  Shock  Zoomed  In 


The  temperature  profile  shows  that  the  SAR  results  seem  to  be  offset  by  about 
18  K  compared  to  Bird  and  NoAR.  The  codes  both  have  the  same  initial  temperature 
in  the  data  subroutine,  so  the  issue  must  be  one  of  sampling. 

An  attempt  to  offset  the  temperature  difference  shown  in  Figure  4.21  by  chang¬ 
ing  the  input  temperature  in  the  data  subroutine  still  does  not  line  up  with  Bird’s 
results.  It  should  be  noted  that  Bird’s  results  matches  the  analytical  temperature 
difference  using  normal  shock  relations,  while  SAR  does  not. 
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Figure  4.21:  Temperature  Line  Plot  for  Mach  1.2  Normal  Shock 


Figure  4.22:  Temperature  Line  Plot  for  Mach  1.2  Normal  Shock  with  Offset 

Additionally,  the  change  in  input  temperature  affects  other  parameters,  as  seen 
in  Figure  4.23.  The  offset  value  causes  the  velocity  after  the  shock  to  be  smaller  than 
it  should  be.  The  issue  with  temperature  needs  to  be  looked  at  further  in  order  to 
understand  why  only  that  flowfield  parameter  is  affected  while  velocity  without  the 
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Figure  4.23:  U-Velocity  Line  Plot  for  Mach  1.2  Normal  Shock  with  Offset 

offset  seems  to  produce  good  results.  A  considerable  amount  of  time  investigating  the 
temperature  issue  has  not  led  to  any  resolution  of  the  problem,  and  DSMC2A  SAR 
algorithms  do  not  have  this  issue. 

Line  plots  have  also  been  produced  for  Mach  9  in  order  to  understand  the  effect 
of  higher  velocities  on  the  flowfield  properties. 

Bird  and  e  =  0%  do  not  match  exactly.  It  was  expected  that  they  would  match 
since  they  have  the  same  accept/reject  criteria  and  density  is  not  sampled  based  on 
weighting.  The  Mach  9  results  show  much  less  variance  than  the  Mach  1.2  results, 
even  though  they  are  both  computed  for  10,000  time  steps,  which  can  be  seen  by  the 
smooth  lines  as  compared  to  Mach  1.2. 
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Figure  4.24:  Density  Line  Plot  for  Mach  9  Normal  Shock 


Figure  4.25:  U-Velocity  Line  Plot  for  Mach  9  Normal  Shock 


Bird  is  not  at  the  extrema,  but  rather  in  the  middle  of  the  SAR  results.  The 
smooth  lines  can  best  be  seen  in  the  magnified  velocity  profile  in  Figure  4.26. 
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Figure  4.26:  U-Velocity  Line  Plot  for  Mach  9  Normal  Shock  Zoomed  In 

Again,  the  SAR  temperature  values  are  offset  from  Bird  and  NoAR  in  Figure 
4.27.  There  does  not  appear  to  be  a  significant  offset  before  the  shock,  but  after  the 
shock  there  is  a  difference  of  about  500  K. 


Figure  4.27:  Temperature  Line  Plot  for  Mach  9  Normal  Shock 
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The  Mach  9  line  plots  are  very  similar  to  Mach  1.2  plots  with  the  exception  of 
the  smoother  lines.  Additionally,  the  change  in  flowheld  properties  across  the  shock 
is  larger,  due  to  the  stronger  shock  typical  of  a  much  higher  Mach  number.  Velocity 
distributions  will  be  evaluated  next  in  order  to  compare  the  macroscopic  properties 
to  the  particle  velocities. 

4-3.2  Velocity  Distributions.  There  are  two  sets  of  velocity  distribution 
data  for  the  SAR  and  NoAR  algorithms.  The  first  set  is  all  the  particles  in  the  cell, 
without  regard  to  weighting,  which  will  show  how  the  surface  would  experience  the 
flow  in  the  SAR  and  NoAR  algorithms,  since  weighting  is  not  taken  into  account  for 
the  reflect  module  [4],  The  weighted  and  unweighted  distributions  will  also  show  if 
the  changes  are  due  more  to  the  change  in  collision  rate  due  to  the  SAR  method  or 
if  they  are  due  to  the  weighting  scheme.  The  velocity  distributions  were  sampled  at 
four  different  Mach  numbers:  1.5,  3,  6,  and  9  before  the  shock,  in  the  shock,  and  after 
the  shock.  The  distributions  will  be  compared  each  other,  and  are  sampled  after  the 
collision  in  order  to  correctly  investigate  the  influence  of  the  accept /reject  changes  on 
the  flowfield. 


Figure  4.28:  Unweighted  Velocity  Distributions  Before  Shock  Mach  1.5 

In  Figure  4.28,  the  unweighted  velocity  distributions  for  Bird,  e  values  of  0%, 
50%  and  100%,  and  NoAR  before  the  shock  was  plotted  in  the  same  manner  as  the 
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simulation  data  compared  to  Holtz  and  Muntz  was  plotted.  The  distributions  are  all 
nearly  identical,  which  shows  again  that  the  change  in  accept/reject  criteria  is  not  the 
main  cause  of  the  flowheld  changes,  but  rather  the  weighting  associated  with  SAR 
created  the  large  differences.  All  of  the  simulation  results  have  a  peak  to  the  right  bulk 
velocity.  Since  these  distributions  are  taken  before  the  shock,  one  would  expect  the 
peak  to  be  at  one.  The  distributions  are  wider  than  the  upstream  velocity  distribution 
in  Figure  2.11,  which  is  due  to  the  lower  velocity  of  the  flow  field.  The  unweighted 
e  =  0%  distribution  is  identical  to  Bird’s  distribution,  as  is  expected,  e  =  0%  has 
the  same  accept/reject  criteria  as  Bird,  therefore  the  only  difference  between  a  SAR 
algorithm  with  e  =  0%  and  DSMC  is  the  weighting.  If  weighting  is  neglected,  the 
flowfields  match  exactly. 


Figure  4.29:  Weighted  Velocity  Distributions  Before  Shock  Mach  1.5 

The  weighted  distributions  in  Figure  4.29  are  jagged  compared  to  the  un¬ 
weighted  distributions,  which  due  smaller  sample  size.  The  peak  of  the  weighted 
distributions  is  in  the  same  location  as  the  unweighted  distribution. 
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Figure  4.30:  Unweighted  Velocity  Distributions  In  Shock  Mach  1.5 

In  Figure  4.30  the  distributions  are  very  similar  to  the  distributions  before  the 
shock.  A  shock  at  Mach  1.5  is  fairly  weak,  and  the  amount  of  nonequilibrium  is  fairly 
small.  Therefore,  the  macroscopic  flowheld  parameters  and  the  velocity  distributions 
will  show  smaller  changes  compared  to  shocks  at  higher  Mach  numbers. 
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Figure  4.31:  Weighted  Velocity  Distributions  In  Shock  Mach  1.5 


In  Figure  4.31,  the  weighted  distributions  are  slightly  less  jagged  than  before 
the  shock,  which  is  due  to  a  larger  number  density  in  the  shock  which  increases  the 
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sample  size  of  the  distribution.  The  peaks  for  the  weighted  distributions  are  in  about 
the  same  location  as  the  peaks  of  the  unweighted  distributions. 


Figure  4.32:  Unweighted  Velocity  Distributions  After  Shock  Mach  1.5 


After  the  shock,  in  Figure  4.32,  the  distributions  again  look  very  similar  to  the 
distributions  before  the  shock  and  in  the  shock.  Since  the  shock  is  weak,  it  is  expected 
that  the  distributions  would  not  change  very  much,  but  it  was  expected  that  the  peak 
of  the  distribution  would  be  at  a  location  less  than  one. 
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Figure  4.33:  Weighted  Velocity  Distributions  After  Shock  Mach  1.5 
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The  weighted  distributions  in  Figure  4.33  are  still  fairly  jagged.  The  changes  in 
the  distribution  across  the  shock  are  very  minor,  so  a  larger  Mach  number  will  need 
to  be  investigated  in  order  to  show  how  the  distributions  compare  when  there  is  a 
larger  gradient  within  the  shock,  which  leads  to  more  of  a  nonequilibrium  condition. 


Figure  4.34:  Unweighted  Velocity  Distributions  Before  Shock  Mach  3 

As  seen  in  Figure  4.34,  at  Mach  3,  the  preshock  distributions  are  thinner  than 
the  Mach  1.5  preshock  distributions.  The  distributions  are  still  centered  to  the  right 
of  the  bulk  velocity,  but  they  are  not  as  far  off  as  they  were  in  the  Mach  1.5  case.  All 
of  the  unweighted  distributions  are  very  similar  to  each  other,  just  as  was  seen  with 
the  Mach  1.5  figures. 
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Figure  4.35:  Weighted  Velocity  Distributions  Before  Shock  Mach  3 


The  weighted  distributions  in  Figure  4.35  are  very  jagged  compared  to  the 
unweighted  distributions.  The  e  =  100%  distribution  is  less  jagged  than  the  e  =  0% 
distribution,  and  NoAR  is  the  least  jagged  of  the  weighted  distributions.  As  e  value 
increases,  so  does  the  sample  size,  and  NoAR  has  the  largest  sample  size  of  weighted 
distributions  and  therefore  is  more  smooth. 


Figure  4.36:  Unweighted  Velocity  Distributions  In  Shock  Mach  3 
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Figure  4.36  shows  a  wider  distribution  for  all  of  the  cases  compared  to  the 
distributions  before  the  shock.  Also,  the  distributions  are  centered  around  a  value 
less  than  the  bulk  velocity,  as  expected. 
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Figure  4.37:  Weighted  Velocity  Distributions  During  Shock  Mach  3 


Just  as  with  Mach  1.5,  the  weighted  distributions  in  the  shock  are  less  jagged 
than  before  the  shock  due  to  an  increase  in  density  in  Figure  4.37.  The  weighted 
NoAR  distribution  is  very  smooth,  and  other  than  being  flatter  at  the  top  of  the 
distribution,  very  closely  resembles  the  unweighted  NoAR  distribution. 


Figure  4.38:  Unweighted  Velocity  Distributions  After  Shock  Mach  3 
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Figure  4.38  shows  that  after  shock  the  distributions  are  centered  to  the  left  of 
the  bulk  velocity,  which  indicates  that  the  particles  have  reequilibrated  around  a  new 
bulk  velocity  that  is  lower  than  the  upstream  bulk  velocity.  So  far,  there  has  been 
very  little  difference  between  the  unweighted  distributions  at  each  location,  which 
indicates  that  the  change  in  the  collision  rate  does  not  play  as  large  of  a  role  in  the 
flow  field  compared  to  the  weighting. 


Figure  4.39:  Weighted  Velocity  Distributions  After  Shock  Mach  3 

The  weighted  distributions  in  Figure  4.39  after  the  shock  are  very  similar  to 
each  other,  as  is  expected  because  this  portion  of  the  flowfield  is  in  equilibrium.  The 
weighted  NoAR  distribution  again  is  the  most  smooth  of  the  weighted  distributions. 
There  is  very  little  difference  between  the  weighted  distributions  in  the  shock  and  the 
distribution  after  the  shock,  except  for  the  fact  that  the  distributions  after  the  shock 
are  more  smooth. 


69 


Figure  4.40:  Unweighted  Velocity  Distributions  Before  Shock  Mach  6 


In  Figure  4.40,  the  unweighted  distributions  are  more  thin  than  the  lower  Mach 
numbers  already  investigated.  These  distributions,  which  are  at  Mach  6,  are  centered 
on  the  bulk  velocity,  unlike  the  lower  Mach  numbers. 
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Figure  4.41:  Weighted  Velocity  Distributions  Before  Shock  Mach  6 


The  weighted  NoAR  is  still  the  most  smooth  of  the  weighted  distributions  in 
Figure  4.41.  The  weighted  e  =  0%  distribution  is  jagged  at  the  top.  The  distributions 
are  centered  on  one,  as  expected. 
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Figure  4.42:  Unweighted  Velocity  Distributions  In  Shock  Mach  6 


In  Figure  4.42  the  unweighted  distributions  have  a  peak  at  the  bulk  velocity, 
but  the  distribution  is  asymmetric  and  skewed  to  the  right.  Particles  have  collided 
and  are  moving  slower  than  the  bulk  velocity,  but  there  are  still  a  significant  portion 
of  particles  traveling  at  the  bulk  velocity. 


Figure  4.43:  Weighted  Velocity  Distributions  In  Shock  Mach  6 


The  weighted  distributions  in  Figure  4.43  has  a  peak  to  the  left  of  the  bulk 
velocity,  which  makes  sense  because  the  sampling  is  only  from  particles  that  have 
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collided.  The  weighted  distributions  are  more  symmetric  and  look  more  like  the 
equilibrium  distributions  seen  outside  of  the  shock  layer,  which  is  expected. 


Figure  4.44:  Unweighted  Velocity  Distributions  After  Shock  Mach  6 


Figure  4.44  shows  that  after  the  shock,  the  unweighted  distributions  again  look 
very  similar  to  each  other,  indicating  that  in  equilibrium  the  algorithms  will  simulate 
similar  flow  fields.  The  peak  of  the  distributions  are  to  the  left  of  one,  which  shows 
that  the  bulk  velocity  after  the  shock  is  less  than  the  bulk  velocity  before  the  shock. 


Figure  4.45:  Weighted  Velocity  Distributions  After  Shock  Mach  6 
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The  weighted  distributions  in  Figure  4.45  are  very  similar  to  the  unweighted 
distributions  shown  in  the  previous  figure.  The  weighted  distributions  are  less  jagged 
than  the  distributions  at  lower  Mach  numbers,  which  is  due  to  the  increased  density 
after  this  strong  shock. 


Figure  4.46:  Unweighted  Velocity  Distributions  Before  Shock  Mach  9 


In  Figure  4.46,  the  distributions  are  all  very  thin  and  identical  to  each  other. 
The  peaks  are  located  at  the  bulk  velocity. 


Figure  4.47:  Weighted  Velocity  Distributions  Before  Shock  Mach  9 
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The  weighted  distributions  in  Figure  4.47  are  thin  and  very  similar  to  the  un¬ 
weighted  distributions  in  Figure  4.46.  NoAR  has  the  smoothest  distribution  compared 
to  the  other  weighted  distributions,  but  overall  they  are  comparable. 


Figure  4.48:  Unweighted  Velocity  Distributions  In  Shock  Mach  9 


The  bimodal  distribution  as  discussed  by  Holtz  and  Muntz  [19]  is  more  distin¬ 
guishable  at  Mach  9  in  the  unweighted  distributions  in  Figure  4.48.  NoAR  shows  the 
least  bimodal  behavior  compared  to  the  other  unweighted  distributions. 


Figure  4.49:  Weighted  Velocity  Distributions  In  Shock  Mach  9 
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The  distributions  in  Figure  4.49  are  not  bimodal  because  all  the  particles  have 
collided  and  are  in  equilibrium  with  each  other.  The  weighted  distributions  match 
each  other  well  and  are  nearly  symmetric  around  the  peak. 


Figure  4.50:  Unweighted  Velocity  Distributions  After  Shock  Mach  9 

After  the  shock,  in  Figure  4.50,  the  unweighted  distributions  return  to  an  equi¬ 
librium  distribution  with  a  peak  at  less  than  the  bulk  velocity.  The  distribution  after 
the  shock  is  much  wider  than  before  the  shock  because  the  particle  collisions  within 
the  shock  tend  to  create  a  distribution  wider  than  before  the  particles  collided. 


Figure  4.51:  Weighted  Velocity  Distributions  After  Shock  Mach  9 
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In  Figure  4.51,  the  weighted  distributions  are  smoother  than  at  lower  Mach 
numbers,  but  still  not  as  smooth  as  the  unweighted  distributions  show  in  the  previous 
figure.  The  distributions  are  very  similar  to  the  weighted  distributions  within  the 
shock,  but  the  peak  is  farther  to  the  left  after  the  shock,  the  nonequilibrium  areas 
are  where  the  major  differences  between  the  algorithms  will  be  observed. 

Based  on  the  Holtz  and  Muntz  comparison,  and  the  further  inspection  of  velocity 
distributions  at  four  different  Mach  numbers  it  appears  that  Bird’s  code  does  not 
allow  for  enough  collision  to  occur  in  order  to  properly  equilibrate  the  cell.  In  order 
to  confirm  this  idea,  the  Mach  9  simulation  using  Bird’s  code  was  run  again,  but 
the  velocities  of  collided  particles  were  accounted  for  separately  from  the  uncollided 
particles.  The  collided  and  uncollided  distributions  were  graphed  both  normalized 
and  not  normalized  in  order  to  allow  the  reader  to  see  the  overall  affect  and  to  be 
able  to  compare  the  distributions  to  distributions  already  shown. 


(a)  Not  Normalized  (b)  Normalized 

Figure  4.52:  Velocity  Profiles  for  Collided  and  Uncollided  Particles  In  Shock  M=9 


As  seen  previously  with  the  Mach  9  distribution  in  the  shock,  there  is  a  definite 
peak  at  the  bulk  velocity.  Looking  at  just  uncollided  particles  results  in  a  distribution 
that  is  very  similar  to  the  distribution  one  gets  when  looking  at  all  of  the  particles, 
regardless  of  weighting.  However,  the  distribution  of  just  the  collided  particles  is 
substantially  different.  The  smaller  peak  is  due  to  particles  that  have  collided  in 
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previous  cells.  DSMC  re-indexes  the  particles  after  each  time  step  because  of  particles 
entering  and  leaving  the  flowfield.  Therefore  it  is  not  currently  possible  to  track 
the  history  of  the  particles  and  make  a  distribution  of  particles  that  have  collided 
previously  and  a  separate  distribution  of  particles  that  have  never  collided.  The 
distribution  that  has  not  been  normalized  shows  that  the  number  of  collided  particles 
is  much  smaller  than  the  number  of  collided  particles.  The  distribution  is  incredibly 
short  compared  to  the  uncollided  distribution,  but  nearly  as  wide.  The  normalized 
distribution  again  shows  the  bimodal  distribution  for  the  uncollided  particles,  but  the 
collided  particles  look  very  similar  to  the  e  =  0%  weighted  distribution  which  has 
been  overlayed  on  the  plot.  Next,  the  statistical  properties  of  the  distributions  looked 
at  so  far  should  be  investigated  for  a  better  comparison. 

The  column  on  the  left  in  Table  4.2  is  the  unweighted  values,  and  the  column  on 
the  right  is  the  weighted  values  with  the  values  grouped  by  algorithm  and  location. 
Kurtosis  and  skewness  are  calculated  using  equations  3.6  and  3.7.  Cmp  is  the  velocity 
that  corresponds  to  the  peak  of  the  distribution,  which  loses  some  meaning  with  the 
bimodal  distributions,  so  the  average  velocity  is  also  calculated.  The  input  velocity 
for  Mach  1.5  is  478  m/s,  and  the  analytically  calculated  value  for  the  velocity  after 
the  shock  is  278  m/s.  As  noticed  with  the  distributions,  the  Cmp  values  for  all  the 
Mach  1.5  distributions  are  higher  than  the  input  velocity.  The  average  velocities  are  a 
little  higher  than  the  most  probable  velocities,  which  implies  that  the  distributions  are 
skewed  to  the  right.  The  skewness  values  are  positive,  which  also  indicates  that  the 
distributions  are  skewed  to  the  right.  The  kurtosis  values  only  vary  slightly  through 
the  shock  as  would  be  expected  given  the  minor  changes  observed  when  visually 
comparing  the  distributions.  The  shock  is  weak  at  Mach  1.5,  therefore  the  flowfield 
does  not  change  significantly  through  the  shock. 
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Table  4.2:  Distribution  Properties  M=1.5 


Bird 

Kurtosis 

Skewness 

Cmp 

QJ 

> 

bn 

:=- 

< 

before 

2.0221 

0.7642 

544.89 

604.23 

during 

1.8668 

0.6671 

518.79 

565.67 

after 

1.9740 

0.7413 

494.94 

547.41 

EpsOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

2.0221 

0.7642 

544.89 

604.23 

during 

1.8668 

0.6671 

518.79 

565.67 

after 

1.9740 

0.7413 

494.94 

547.41 

Eps25 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

1.9627 

0.7266 

574.44 

604.24 

0.6391 

528.21 

566.72 

0.8558 

500.10 

547.69 

Eps50 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

1.8073 

0.6025 

589.31 

604.57 

during 

1.8174 

0.6308 

522.26 

566.27 

after 

1.9542 

0.7266 

456.59 

547.81 

Eps75 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

1.8271 

0.6288 

572.24 

604.08 

during 

1.7691 

0.5896 

532.97 

566.37 

after 

1.9076 

0.6934 

500.12 

547.17 

EpslOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

1.8643 

0.6445 

586.36 

604.23 

during 

1.8222 

0.6318 

514.14 

566.08 

after 

1.9807 

0.7446 

475.50 

547.04 

NoAR 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

1.9125 

0.6859 

570.19 

603.80 

during 

1.8467 

0.6516 

556.45 

565.85 

after 

2.0548 

0.7990 

508.00 

547.25 

(a)  Unweighted 


EpsOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

2.0394 

0.7419 

604.37 

617.08 

during 

1.8794 

0.6571 

512.27 

578.63 

after 

1.9532 

0.7213 

460.89 

562.21 

Eps25 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

2.0647 

0.7295 

568.00 

611.97 

during 

1.7927 

0.5937 

560.34 

577.21 

after 

2.1556 

0.8416 

537.07 

563.01 

Eps50 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

1.9954 

0.6500 

512.63 

616.19 

during 

1.7965 

0.6036 

490.44 

573.87 

after 

1.9090 

0.6933 

449.86 

559.62 

Eps75 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

1.9166 

0.6429 

614.57 

608.73 

during 

1.7417 

0.5492 

551.64 

576.95 

after 

1.9057 

0.6714 

480.42 

561.00 

EpslOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

1.7874 

0.5720 

519.87 

612.48 

during 

1.8353 

0.6169 

571.38 

575.84 

after 

1.9693 

0.7147 

509.53 

556.88 

NoAR 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

1.9674 

0.7016 

614.02 

601.82 

during 

1.9672 

0.7009 

550.01 

565.30 

after 

2.0585 

0.7889 

458.41 

546.61 

(b)  Weighted 


At  Mach  3,  The  input  bulk  velocity  is  956  m/s,  which  matches  best  with  the 
Cmp  for  the  weighted  e  =  75%.  The  unweighted  Cmp  values  are  all  above  the  input 
velocity,  as  are  the  average  velocity  values,  e  =  0%  and  75%  have  Cmp  values  below 
the  input  velocity,  but  the  rest  of  the  weighted  values  are  above  the  input  velocity. 
Kurtosis  values  are  above  the  kurtosis  values  at  Mach  1.5,  which  is  expected  because 
the  distributions  become  thinner  as  Mach  number  increases.  The  skewness  values 
before  the  shock  are  higher  than  Mach  1.5  skewness  values,  indicating  that  there  are 
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Table  4.3:  Distribution  Properties  M=3 


Bird 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

2.2735 

0.9061 

1012.04 

1020.24 

during 

1.7173 

0.5669 

336.47 

367.43 

after 

2.0991 

0.8183 

717.74 

310.33 

EpsOO 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.2735 

0.9061 

1012.04 

1020.24 

during 

1.7173 

0.5669 

336.47 

367.43 

after 

2.0991 

0.8183 

717.74 

310.33 

Eps25 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.2041 

0.8693 

993.99 

1020.01 

during 

1.7939 

0.6283 

902.10 

367.13 

after 

2.0020 

0.7573 

713.47 

310.05 

Eps50 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.1292 

0.8214 

1021.60 

1020.61 

during 

1.7081 

0.5576 

303.81 

367.44 

after 

1.9053 

0.6917 

723.40 

310.53 

Eps75 

Kurtosis 

Skewness 

before 

2.1734 

0.8536 

1015.54 

1020.06 

during 

1.7374 

0.5333 

362.12 

363.19 

after 

2.0173 

0.7633 

737.52 

310.34 

EpslOO 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

n 

0.8391 

1059.76 

1020.34 

during 

1.6658 

0.5155 

337.15 

363.35 

after 

1.9002 

0.6893 

743.11 

309.97 

NoAR 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.2751 

0.9130 

935.30 

1020.60 

during 

1.9210 

0.7261 

366.13 

363.60 

after 

1.9625 

0.7324 

713.15 

309.64 

EpsOO 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.2935 

0.3694 

973.34 

1028.29 

during 

1.6937 

0.5313 

321.83 

330.43 

after 

2.0575 

0.7317 

664.43 

334.05 

Eps25 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.1639 

0.3147 

946.06 

1022.38 

during 

1.7314 

0.6013 

777.59 

376.20 

0.7209 

749.10 

332.71 

Eps50 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.0127 

0.7399 

970.75 

1021.79 

during 

1.6356 

0.5314 

317.97 

330.94 

after 

1.3425 

0.6473 

694.03 

334.73 

Eps75 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.3998 

0.9179 

943.15 

1025.41 

during 

1.7234 

0.5603 

390.15 

331.62 

after 

1.9668 

0.7309 

686.00 

334.12 

EpslOO 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.2330 

0.3314 

1007.92 

1024.84 

during 

1.7015 

0.5193 

823.22 

373.93 

after 

1.3340 

0.6443 

733.34 

332.67 

NoAR 

Kurtosis 

Skewness 

Cmp 

AvgVel 

before 

2.2931 

0.9259 

977.51 

1022.31 

during 

1.9245 

0.7203 

753.20 

366.51 

after 

1.9591 

0.7320 

723.20 

310.37 

(a)  Unweighted  (b)  Weighted 

more  particles  traveling  faster  than  Cmp  at  Mach  3  than  Mach  1.5.  The  skewness  of 
the  distribution  in  the  shock  decreases  compared  to  before  the  shock,  and  increases 
again  after  the  shock. 

In  Table  4.4,  the  kurtosis  is  higher  than  the  lower  Mach  numbers,  as  is  expected. 
The  skewness  for  all  the  results  is  very  high  before  the  shock,  and  within  the  shock 
drops  to  about  half  the  before  shock  value.  The  unweighted  distributions  are  begin¬ 
ning  to  exhibit  more  of  a  bimodal  distribution  within  the  shock,  with  the  majority  of 
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Table  4.4:  Distribution  Properties  M=6 


Bird 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

3.0933 

1.2717 

1946.33 

1944.65 

during 

1.3933 

0.7060 

1315.79 

1583.58 

after 

2.0761 

0.3065 

1238.35 

1453.75 

before 


Avg  ve 


1944.65 


3 

0.7060 

1315.79 

1583.5 

1288.85  1453.75 


Eps25 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

2.6076 

1.0729 

1953.33 

1944.69 

during 

1.3365 

0.7025 

1307.26 

1533.67 

after 

1.9530 

0.7233 

1254.33 

1454.14 

Eps50 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

2.7235 

1.1255 

1942.74 

1945.64 

during 

1.7612 

0.6013 

1332.17 

1583.42 

after 

2.1515 

0.3533 

1214.43 

1453.14 

Eps75 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

4.3638 

1.6362 

1925.66 

1944.99 

during 

1.3178 

0.6536 

1319.05 

1590.00 

after 

2.0193 

0.7705 

1293.35 

1453.81 

EpslOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

3.0760 

1.2700 

1957.33 

1944.77 

during 

1.3736 

0.6971 

1357.43 

1592.33 

after 

1.9923 

0.7532 

1241.34 

1453.60 

NoAR 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

3.1639 

1.2992 

1957.49 

1945.07 

during 

1.3043 

0.6449 

1735.66 

1592.13 

after 

2.0344 

0.3123 

1239.33 

1452.84 

EpsOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

2.8997 

1.2031 

1395.83 

1942.54 

during 

1.8472 

0.6502 

1431.42 

1597.09 

after 

2.0121 

0.7634 

1307.30 

Em 

Eps25  kurtosis  Skewness 


before  2.5004  1.0153 

during  1.8285  0.6436 

after  1.8882  0.6871 


;ewness  Ctnp  AvgV 


1.0153  1360.66  1947. 

0.6436  1578.31  1606. 

0.6871  1344.94  1497. 


Eps50  Kurtosis  Skewness  Cmp 

before  2.8065  1.1458  1396.43 

during  1.6303  0.5348  1503.73 

after  2.0969  0.3212  1352.09 


1943.76 

1603.75 


Eps75 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

4.3131 

1.6639 

1992.65 

1943.26 

during 

1.7670 

0.6030 

1562.42 

1607.37 

after 

1.9663 

0.7336 

1238.21 

1496.06 

100  Kurtosis  Skewness 
ore  3.0462  1.2378 

■ing  1.8343  0.6540 


after  1.9450  0.7225 


Cmp 

1968.04 

1594.47 


1369.47 


Avg  Vel 
1949.94 
1602.79 


1496.01 


NoAR  kurtosis  S 


before  3.3080 
during  1.7659  i 


1.3231  1926.31  1947.13 
0.6102  1718.66  1584.39 


Unweighted 


0.3141  1270.79  1451.10 


(b)  Weighted 


particles  at  the  bulk  velocity,  which  results  in  a  slightly  higher  skewness  compared  to 
the  weighted  distributions.  The  Cmp  values  are  closer  to  the  bulk  velocity  compared 
to  lower  Mach  numbers  with  the  unweighted  e  =  75%  value  being  closest  to  1913 
m/s  followed  closely  by  weighted  NoAR.  The  average  velocities  match  closely  to  the 
corresponding  Cmp  for  each  case,  which  is  most  likely  due  to  the  thin  distributions  at 
this  Mach  number. 
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Table  4.5:  Distribution  Properties  M=9 


Bird 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

7.9963 

2.5040 

2920.43 

2889.46 

during 

2.2365 

0.8832 

2836.36 

2342.86 

after 

1.9950 

0.7557 

1879.09 

2132.22 

EpsOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

8.0031 

2.5052 

2920.43 

2889.90 

during 

2.2259 

0.8790 

2840.32 

2342.91 

after 

2.0709 

0.8040 

1885.12 

2131.77 

Eps25 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

8.0031 

2.5052 

2920.43 

during 

2.2259 

0.8790 

2840.32 

2342.91 

after 

2.0709 

0.8040 

1885.12 

2131.77 

Eps50 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

5.4030 

1.9559 

2861.51 

2891.05 

during 

1.9425 

0.7029 

2823.89 

2342.08 

after 

2.2519 

0.9091 

1789.82 

2131.85 

Eps75 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

5.6947 

2.0239 

2922.12 

2890.30 

during 

2.0237 

0.7615 

2840.26 

2344.60 

after 

1.9517 

0.7242 

1830.55 

2129.80 

EpslOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

5.6799 

2.0230 

2892.93 

2890.86 

BBS 

2731.01 

2345.75 

after 

2.1727 

0.8672 

1971.78 

2133.05 

NoAR 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

5.9486 

2.0822 

2895.75 

2890.21 

during 

1.8620 

0.6631 

2737.73 

2349.70 

after 

2.0118 

0.7633 

1873.70 

2131.64 

(a)  Unweighted 


EpsOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

8.0371 

2.5070 

2854.28 

2894.51 

during 

1.9988 

0.7696 

2312.20 

2362.23 

after 

2.0230 

0.7711 

1857.28 

2193.72 

Eps25 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

8.0371 

2.5070 

2894.51 

during 

1.9988 

0.7696 

2362.23 

after 

2.0230 

0.7711 

2193.72 

Eps50 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

5.3606 

1.9472 

2813.49 

2893.92 

during 

1.7342 

0.5831 

2453.08 

2359.71 

after 

2.2049 

0.8792 

1999.04 

2193.51 

Eps75 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

before 

5.3905 

1.9572 

2938.74 

2892.81 

during 

1.8209 

0.6484 

2247.58 

2359.73 

after 

1.9092 

0.6946 

1961.83 

2190.56 

EpslOO 

Kurtosis 

Skewness 

Avg  Vel 

before 

5.4770 

1.9643 

2909.62 

bees 

0.7133 

2486.22 

2359.85 

after 

2.1322 

0.8413 

1942.71 

2194.89 

NoAR 

Kurtosis 

before 

5.9833 

2.0894 

2893.26 

during 

1.7779 

0.6049 

2323.30 

2346.14 

after 

2.0070 

1846.73 

2133.53 

(b)  Weighted 


The  before  shock  values  for  kurtosis  and  skewness  are  much  higher  compared 
to  previous  Mach  numbers  and  the  weighted  and  corresponding  unweighted  distribu¬ 
tions  have  very  similar  values.  Within  the  shock,  skewness  is  again  higher  for  the 
unweighted  distributions.  The  Cmp  values  before  the  shock  are  fairly  close  to  the  av¬ 
erage  velocities,  with  the  exception  of  the  SAR  cases  previously  discussed.  The  Cmp 
value  closest  to  the  input  velocity  of  2867  m/s  is  the  weighted  NoAR  value  of  2878 
m/s.  The  bimodal  distribution  seen  with  the  unweighted  distributions  leads  to  a  lower 
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average  velocity  compared  to  the  most  probable  velocity.  Evaluating  these  velocity 
distributions  have  allowed  for  a  better  understanding  of  the  microscopic  behavior  of 
the  particles,  and  the  change  of  the  behavior  when  using  SAR  and  NoAR.  Through 
the  shock,  the  number  of  collisions  as  calculated  in  Bird’s  code  is  obviously  too  low, 
which  causes  a  longer  equilibration  thereby  making  the  shock  layer  thicker  compared 
to  experimental  data. 
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4-4  2d  Axisymmetric  Hollow  Cylinder  Results 

The  DSMC2A  code  was  used  for  a  variety  of  analyses:  density  plots  compared 
to  Davis’  experimental  data,  percent  difference  contour  plots,  surface  plots,  contour 
plots,  and  the  cell  size  was  reduced  to  one  quarter  the  original  size  with  results  were 
compared  to  Davis’  data.  Before  viewing  the  results  in  comparison  to  experimental 
data,  the  contour  plots  should  be  investigated  so  that  the  overall  flowheld  is  under¬ 
stood.  Since  density  is  the  primary  flowfield  parameter  investigated  in  this  project, 
it  should  be  the  first  one  discussed.  There  are  only  slight  variations  in  the  flowfields 
between  Bird,  SAR,  and  NoAR,  so  only  Bird,  e  =  0%,  e  =  100%,  and  NoAR  will  be 
shown  in  this  section. 
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Figure  4.53:  Density  Contour  Plots 


The  stagnation  point  is  very  visible  in  Figure  4.53  for  all  the  cases,  as  is  the 
shock  layer.  Only  very  slight  variations  in  the  shock  layer  can  be  seen  when  comparing 
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the  cases.  NoAR  results  in  a  shock  layer  closer  to  the  step  that  appears  to  be  thicker 
than  Bird’s  shock  layer.  The  temperature  plots  are  also  very  similar  to  each  other  in 
Figure  4.54. 


(c)  EpslOO  (d)  NoAR 


Figure  4.54:  Translational  Temperature  Contour  Plots 

An  area  of  hot  gas  can  be  seen  before  the  step  in  Figure  4.54,  where  the  shock 
layer  and  boundary  layer  are  converging.  The  size  and  temperature  of  that  area  varies 
with  algorithm.  The  shock  layer  above  the  step  can  also  be  seen,  with  minor  changes 
for  each  case.  Note  that  on  the  horizontal  face  of  the  step  in  Bird’s  case,  the  step  does 
not  seem  flat:  there  are  pockets  of  hot  air  throughout  the  length  of  the  step.  The 
issue  on  the  step  is  due  to  low  populations  of  particles  in  the  cells  in  that  area,  which 
causes  statistical  errors.  The  SAR  and  NoAR  cases  do  not  suffer  the  same  problem, 
even  though  they  also  have  few  cells  behind  the  shock  along  the  horizontal  face  of  the 
step,  which  will  be  shown  shortly. 
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(c)  EpslOO 


(d)  NoAR 


Figure  4.55:  U-Velocity  Contour  Plots 


In  Figure  4.55a,  the  same  problem  along  the  top  of  the  step  can  be  seen  again, 
but  it  is  not  seen  with  SAR  and  NoAR.  In  order  to  show  the  amount  of  rarefaction  in 
each  cell,  these  cases  were  run  again  and  the  number  of  particles  per  cell  was  printed 
out  to  a  hie. 

As  seen  in  Figure  4.56,  the  number  of  particles  is  less  than  10  on  the  horizontal 
face  of  the  step.  In  the  stagnation  region,  the  number  rises  to  about  40  particles,  and 
in  the  shock  layer  the  number  of  particles  increases  even  more  to  approximately  60 
particles.  The  howheld  makes  choosing  a  ratio  of  simulated  particles  to  real  particles  a 
difficult  one.  If  the  ratio  is  too  small,  the  rarefied  areas  of  the  how  become  even  more 
rarehed,  leading  to  errors.  If  the  ratio  is  too  large,  particles  will  collide  more  than 
once  per  time  step  in  stagnation  regions  and  in  the  shock  layer  and  the  assumption  of 
decoupled  particle  collisions  and  motion  is  no  longer  valid.  SAR  and  NoAR  provide 
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Figure  4.56:  Average  Number  of  Particles  Per  Cell 


good  results,  even  in  areas  of  higher  rarefaction,  which  is  an  important  point  to 
consider.  The  sampling  algorithm  in  SAR  and  NoAR  is  most  likely  what  causes 
the  change,  since  e  =  0%  results  do  not  have  the  same  problem  as  Bird’s  results  even 
though  their  accept/reject  criteria  are  identical.  The  weighting  algorithm  in  SAR  and 
NoAR  most  likely  changes  the  variance  significantly  in  these  low  populated  regions, 
and  allows  for  the  flowfield  properties  to  be  correctly  sampled. 
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4-4-1  Comparison  to  Davis’  Experimental  Data.  Previous  comparison  to 
experimental  data  showed  that  the  NoAR  and  Bird  results  were  usually  the  extrema 
with  the  SAR  values  in  between  with  the  0%  case  closest  to  Bird  and  higher  e  values 
tending  toward  NoAR  [11].  The  density  profile  results  from  the  current  simulation 
shows  the  same  trend. 


Figure  4.57:  Density  Profile  at  x=. 0313m 


In  Figure  4.57,  the  profile  data  is  taken  .0187m  upstream  of  the  step.  All  the 
cases  result  in  a  shape  that  is  fairly  similar  to  the  experimental  data,  but  SAR  with 
e  =  50%  is  matches  best  with  the  change  in  density.  Throughout  the  profile  Bird  and 
e  =  00%  match,  which  is  why  Bird’s  green  line  does  not  appear. 
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Figure  4.58:  Density  Profile  at  x=. 0462m 

Again,  all  the  cases  provide  results  that  have  the  same  shape,  but  e  =  50%  has 
the  same  change  in  density  at  y=0.032  m.  None  of  the  computational  results  match 
the  solution  closer  to  the  wall,  but  NoAR  and  the  high  e  value  of  200%  are  the  closest. 


Figure  4.59:  Density  Profile  at  x=. 0495m 


At  the  step,  the  e  =  25%  value  is  closest  to  the  density  peak  at  y=  0.033  m, 
and  NoAR  and  e  =  200%  match  the  best  closest  to  the  wall.  All  the  computational 
results  match  the  overall  behavior  of  the  density  profile. 


Figure  4.60:  Density  Profile  at  x=. 0509m 

Past  the  stagnation  point  in  Figure  4.60,  the  density  profile  is  similar  to  the 
first  profile  at  x=0.0313  m.  The  change  in  density  is  larger,  as  the  shock  is  stronger 
after  the  step.  Again,  just  after  the  step  e  =  25%  matches  best  with  the  change  in 
density. 
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Figure  4.61:  Density  Profile  at  x=0. 0561m 

The  shock  is  even  stronger  6  mm  downstream  from  the  step,  as  shown  by  the 
change  in  density  at  x=0.0561  m  in  Figure  4.61.  Farther  downstream  from  the  step, 
the  change  in  density  most  closely  matches  with  e  =  100%.  Below  y=0.034  m,  Bird 
and  e  =  0%  match  the  best  to  the  experimental  plots.  None  of  the  plots  show  good 
agreement  near  the  wall,  which  may  be  due  to  the  specular  boundary  condition  in 
DSMC2A.  At  a  molecular  level,  the  surface  of  the  cylinder  will  not  be  completely  flat, 
which  means  that  a  specular  boundary  condition  is  not  realistic. 

In  an  effort  to  see  why  different  e  values  match  at  different  points  on  the  profile, 
temperature  and  velocity  plots  have  been  created.  These  plots  can  be  compared  to 
the  density  plots  for  an  overall  understanding  of  the  fluid  at  the  sampling  locations. 
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(a)  Velocity 


(b)  Temperature 


Figure  4.62:  Velocity  and  Temperature  Profiles  at  x=0. 0313m 


The  shock  layer  is  not  very  visible  in  the  velocity  profile,  which  matches  the 
velocity  contour  plot.  The  boundary  layer  is  very  clear  in  the  velocity  profile,  however. 
The  temperature  profile  shows  a  slight  change  in  the  temperature  gradient  at  0.029 
m,  which  is  the  peak  of  the  density  profile.  The  temperature  gradient  changes  again 
at  approximately  y=0.025  m,  which  is  where  an  inflection  point  in  the  density  profile 
exists. 


(a)  Velocity  (b)  Temperature 

Figure  4.63:  Velocity  and  Temperature  Profiles  at  x=0. 0462m 
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The  velocity  profile  at  x=0.0462  m  has  the  beginnings  of  what  looks  to  be  an 
unfavorable  gradient,  which  occurs  before  the  boundary  layer  separates  from  the  wall, 
or  in  the  case  of  this  experiment,  as  the  flow  moves  toward  a  stagnation  point.  The 
temperature  profile  has  an  inflection  point  at  y=0.03  m,  which  matches  one  of  the 
inflection  points  in  the  density  profile. 


(a)  Velocity  (b)  Temperature 

Figure  4.64:  Velocity  and  Temperature  Profiles  at  x=. 0495m 

The  velocity  profile  very  clearly  shows  a  stagnation  region  below  y=0.028  m. 
Again,  the  inflection  point  in  the  temperature  profile  corresponds  to  the  second  in¬ 
flection  point  in  Figure  4.59. 
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(a)  Velocity 


(b)  Temperature 


Figure  4.65:  Velocity  and  Temperature  Profiles  at  x=. 0509m 


After  the  step,  the  velocity  profile  returns  to  a  normal  boundary  layer  velocity 
profile.  Interestingly,  the  temperature  profile  in  Figure  4.65  has  an  inflection  point 
that  does  not  match  the  simulated  results,  but  it  does  match  the  experimental  inflec¬ 
tion  point. 


(a)  Velocity 


(b)  Temperature 


Figure  4.66:  Velocity  and  Temperature  Profiles  at  x=. 0561m 
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The  last  velocity  profile  shows  a  higher  boundary  layer,  and  the  temperature 
profile  has  an  inflection  point  that  again  matches  the  density  profile’s  inflection  point. 

All  of  the  plots  show  that  the  computational  results  match  the  same  behavior 
as  the  experimental  profiles.  The  computational  results  tend  to  be  approximately 
0.003  m  higher  than  the  experimental  results,  which  could  easily  be  explained  by 
uncertainty  of  the  experimental  data.  Uncertainty  of  the  experimental  data  and 
the  conversion  of  the  graphical  data  from  Davis’  thesis  to  plots  is  significant,  and 
therefore  the  comparison  of  the  experimental  data  to  the  simulations  is  considered 
reasonable.  However,  due  to  the  uncertainty  associated  with  the  experimental  data,  it 
it  very  difficult  to  definitively  say  that  one  algorithm  is  better  than  the  other.  Rather, 
the  DSMC2A  results  should  be  used  to  show  that  the  SAR  and  NoAR  algorithms 
display  similar  behavior  in  a  2-dimensional  axisymmetric  program  as  they  do  in  the 
1-dimensional  shock  program. 

4-4-2  Quarter  Cell  Size  Cases.  In  Bird’s  book,  it  is  suggested  that  the  cell 
size  be  on  the  order  of  the  mean  free  path  [4].  The  cell  size  for  the  original  case  is 
approximately  1  mm  by  1  mm  while  the  mean  free  path  is  0.47  mm.  The  cell  sized 
was  reduced  to  0.5  mm  by  0.5  mm  to  see  if  the  results  compared  to  Davis  can  be 
further  refined. 
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Figure  4.67:  Density  Profile  at  x=. 0313m 

The  first  profile  shows  that  the  SAR  value  that  would  compare  best  to  the 
experimental  data  is  between  e  values  of  25%  and  50%,  which  is  slightly  lower  than 
the  original  results. 


Figure  4.68:  Density  Profile  at  x=. 0462m 
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Figure  4.68  compares  best  with  25%  at  the  peak  which  is  again  less  than  the 
original  grid.  NoAR  is  closest  below  the  peak. 


0.0001 


□ 

□ 

ft 

,  D^l 


0.0002  0.0003  0.0004 

Density,  kg/mA3 


Davis 

Bird 

EpsOO 

Eps25 

Eps50 

Eps75 

EpslOO 

Eps200 

NoAR 


j  i  i  i  I  i  L 


Figure  4.69:  Density  Profile  at  x=. 0495m 


Figure  4.69  compares  best  with  Bird  and  00%  at  the  peak,  and  below  the  peak 
NoAR  compares  the  best.  Note  that  the  exaggerated  behavior  of  the  computational 
results  compared  to  experimental  results  still  exists. 
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Figure  4.70:  Density  Profile  at  x=. 0509m 


Again,  Bird  and  00%  match  best  with  the  profile  just  after  the  step,  while  25% 
matches  best  for  the  original  grid.  Below  the  peak,  more  of  the  cases  are  closer  to 
experimental  data  compared  to  Figure  4.60. 
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Figure  4.71:  Density  Profile  at  x=. 0561m 
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The  last  profile,  which  is  taken  well  downstream  of  the  step,  matches  closest 
with  50%  at  the  peak,  and  NoAR  closer  to  the  wall.  The  original  grid  matches  closest 
with  e  =  100%. 

The  smaller  grid  consistently  required  a  smaller  e  value  to  match  the  results. 
Since  the  ratio  of  simulated  particles  to  real  particles  was  increased  proportionally  to 
the  decreased  size  of  the  cells,  the  collision  ratio  and  variance  stayed  approximately 
the  same.  This  study  consists  of  comparing  to  experimental  data  with  one  Mach 
number,  which  does  not  allow  for  further  research  into  the  Mach  number  dependency 
for  the  e  input  value.  The  important  lesson  from  this  study  is  that  reducing  the  cell 
size  to  the  length  of  the  mean  free  path  does  significantly  change  the  e  value  that 
matches  experimental  data  through  the  shock  layer.  The  addition  of  more  cells  and 
increase  in  the  ratio  of  particles  caused  the  computational  time  to  grow  significantly. 
Each  of  these  cases  took  approximately  3  days  to  complete,  compared  to  less  than 
half  a  day  for  the  original  grid.  The  smaller  grid  will  provide  a  more  refined  solution, 
but  may  not  be  worth  the  additional  computational  time  to  get  the  better  solution. 

4-4-3  Percent  Difference  Contour  Plots.  Contour  plots  of  the  flowfields 
have  been  shown  earlier  in  this  section  in  order  to  show  the  differences  between  the 
results  of  the  cases.  It  is  difficult  to  see  all  the  differences  with  these  contour  plots, 
however.  The  percent  difference  plots  allow  for  a  better  understanding  of  the  changes 
in  the  flowfield  for  SAR  and  NoAR  compared  to  Bird’s  output.  The  percent  difference 
contour  plots  were  created  for  density,  temperature,  and  u-velocity.  The  contour  plots 
of  percent  difference  for  density  for  each  of  the  cases  is  plotted  in  Figure  4.72. 

In  Figure  4.72a,  all  the  cells  have  a  value  of  zero,  indicating  that  the  flowfield 
exactly  matches  that  of  Bird,  which  make  sense  since  density  is  not  sampled  using  the 
weightings  from  SAR  and  they  have  the  same  accept/reject  criteria.  There  are  large 
differences  in  density  in  the  shock  layer,  and  along  the  horizonontal  face  of  the  step, 
where  issues  with  Bird’s  results  have  previously  been  noted.  The  differences  increase 
with  higher  e  values,  with  NoAR  providing  the  largest  differences.  The  areas  of  large 
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Figure  4.72:  Percent  Difference  Density  Contour  Plots  Compared  to  Bird 


differences  along  the  horizontal  face  of  the  step  is  due  to  an  issue  with  DSMC.  The 
number  of  particles  in  the  area  behind  the  shock  along  the  step  is  very  small,  which 
results  in  a  large  amount  of  variance  and  causes  problems  with  the  correct  simulation 
at  that  location.  Temperature  is  investigated  next. 

The  largest  temperature  differences  in  Figure  4.73  are  found  at  the  leading  edge 
of  the  hollow  cylinder  where  the  shock  layer  begins  and  throughout  the  shock  layer. 
As  the  e  value  increases,  so  do  the  differences  between  Bird  and  the  SAR  case.  NoAR 
shows  the  greatest  difference  compared  to  Bird.  The  SAR  results  show  negative 
differences,  especially  near  the  leading  edge  of  the  cylinder,  while  the  NoAR  results 
does  now  show  any  negative  differences.  The  red  rectangle  is  the  step  that  is  attached 
to  the  hollow  cylinder.  Note  there  is  not  a  large  difference  in  the  boundary  layer, 
even  though  it  is  a  region  of  nonequilibrium.  The  amount  of  nonequilibrium  in  the 
boundary  layer  is  much  smaller  than  the  amount  found  in  the  shock  associated  with 
this  ffowfield. 
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Figure  4.73:  Percent  Difference  Temperature  Contour  Plots  Compared  to  Bird 


Most  of  the  flowfield  matches  Bird’s  results  well  for  u-velocity,  and  the  largest 
differences  are  found  at  the  stagnation  point.  Figure  4.75  focuses  on  the  stagnation 
region  in  order  to  better  view  the  differences  in  that  area. 

Interestingly,  NoAR  shows  the  least  amount  of  difference  compared  to  Bird  in 
the  stagnation  region  in  Figure  4.75.  The  SAR  values  all  show  very  similar  differences, 
but  e  =  0%  shows  the  most  difference  and  e  =  100%  shows  the  least  of  the  SAR  cases. 
The  two  changes  made  to  the  code  involve  the  accept/reject  criteria,  which  affects 
the  collision  rate,  and  the  flowfield  sampling.  Flowfield  sampling  does  not  affect  the 
surface  values,  so  the  percent  difference  for  surface  properties  should  also  be  looked 
at  to  evaluate  the  changes  to  the  system  due  to  the  accept/reject  criteria  alone. 

As  would  be  expected,  in  Figure  4.76,  the  e  =  0%  case  shows  no  difference 
from  Bird’s  results.  The  only  difference  between  Bird  and  e  =  0%  is  the  sampling, 
but  since  sampling  doesn’t  affect  the  surface  values,  they  have  the  same  results.  Both 
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Figure  4.74:  Percent  Difference  U-Velocity  Contour  Plots  Compared  to  Bird 


plots  show  that  NoAR  is  the  most  different  from  Bird  and  the  SAR  results  are  between 
Bird  and  NoAR.  The  vertical  face  shows  a  maximum  difference  of  about  16%  and  the 
horizontal  face  shows  a  maximum  difference  of  18%.  The  differences  on  the  surface 
are  much  smaller  than  the  difference  in  the  temperature  flowfield  values,  which  range 
from  -160%  to  80%,  which  proves  that  sampling  has  a  more  more  profound  effect  on 
the  results  of  the  simulation  than  SAR  does. 


4-4-4  Surface  Plots.  Surface  plots  for  incident  pressure,  incident  transla¬ 
tional  temperature,  and  heat  flux  have  also  been  graphed  for  both  the  horizontal  and 
vertical  faces  of  the  step.  These  plots  can  be  used  to  help  understand  the  effect  of 
the  SAR  and  NoAR  algorithms  on  the  surfaces. 
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(c)  EpslOO  (d)  NoAR 

Figure  4.75:  Percent  Difference  U- Velocity  Contour  Plots  Compared  to  Bird 


The  vertical  face  of  the  step  causes  a  stagnation  point  where  it  meets  the  edge 
of  the  hollow  cylinder.  Heat  flux  is  a  function  of  the  temperature  gradient,  so  it  makes 
sense  that  the  heat  flux  in  Figure  4.77c  follows  the  profile  of  the  incident  temperature 
in  Figure  4.77b. 
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(a)  Vertical  Step  Surface  (b)  Horizontal  Step  Surface 

Figure  4.76:  Surface  Percent  Difference  Plots 


The  horizontal  surface  shows  a  slight  increase  in  incident  pressure  just  after 
the  step,  followed  by  a  steady  decrease  until  x=0.08m  when  there  is  a  sudden  spike 
followed  by  drop  in  pressure.  The  same  behavior  can  be  seen  in  all  of  the  profiles  in 
Figure  4.78. 


4-4-5  Velocity  Distributions.  Now  that  the  macroscopic  properties  have 
been  discussed,  the  velocity  distributions  for  Mach  numbers  1,  3,  6,  and  9  will  be 
shown.  The  DSMC2A  cases  cannot  be  compared  to  a  theoretical  distribution  like 
the  Id  cases.  The  2-dimensional  axisymmetric  flowheld  is  much  more  complicated, 
with  a  shock  layer,  a  circulation  region,  a  stagnation  region,  and  a  boundary  layer. 
These  flowheld  effects  are  associated  with  a  change  in  bulk  velocity,  and  there  is  not 
a  way  to  calculate  what  the  bulk  velocity  should  be  in  each  of  these  cells.  Therefore, 
the  SAR  and  NoAR  speed  distributions  will  be  compared  to  Bird  only,  without  the 
theoretical  distribution  at  each  of  the  sampled  cells.  The  sample  cells  are  located  in 
the  stagnation  region,  boundary  layer,  shock  layer,  and  behind  the  shock. 

As  with  the  1-dimensional  results,  the  unweighted  distributions  are  on  the  left 
and  the  weighted  are  on  the  right  for  all  of  the  figures.  Bird’s  distribution  is  plotted 
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Figure  4.77:  Surface  Plot  for  Front  Face  of  Step 


in  the  background  of  all  of  the  plots  for  comparison.  In  Figures  4.79  and  4.80,  the 
distributions  were  taken  at  the  stagnation  region.  The  stagnation  region  show  particle 
velocities  that  are  slower  than  the  bulk  velocity,  which  is  2000  m/s,  so  the  distribution 
is  to  the  left  of  the  bulk  velocity.  The  unweighted  e  =  0%  distribution  is  nearly 
identical  to  Bird,  which  is  expected  since  they  have  the  same  accept/reject  criteria. 
The  weighted  e  =  0%  distribution  shows  only  slight  differences  compared  to  Bird,  as 
do  the  weighted  and  unweighted  e  =  50%  distributions. 
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Figure  4.78:  Surface  Plot  for  Horizontal  Face  of  Step 

In  Figure  4.80,  the  e  =  100%  distributions  are  again  very  similar  to  Bird,  and 
the  NoAR  distributions  are  the  most  different  from  Bird,  as  would  be  expected. 
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(a)  Unweighted  EpsOO 


(b)  Weighted  EpsOO 


Unweighted  Particle  Velocity  at  N=350  -Eps50 


Weighted  Particle  Velocity  at  N=350  -Eps50 


(c)  Unweighted  Eps50 


(d)  Weighted  Eps50 


Figure  4.79:  Velocity  Distributions  for  DSMC2A  Case  in  Stagnation  Region 


Figure  4.81  shows  a  slightly  wider  distribution  than  was  seen  in  the  stagnation 
region,  and  the  weighted  e  =  0%  and  both  of  the  e  =  50%  cases  are  fairly  similar 
to  Bird.  All  of  the  distributions  are  more  jagged  than  what  was  observed  with  the 
1-dimensional  cases,  indicating  that  the  number  of  particles  are  fewer  in  this  sample 
cell.  The  peaks  of  the  distributions  before  the  step  appear  to  be  farther  to  the  right 
compared  to  the  distributions  at  the  stagnation  point,  which  means  that  the  particles 
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(a)  Unweighted  EpslOO 


(b)  Weighted  EpslOO 


Onweighted  Particle  Velocity  at  N=350  -NoAR  Weighted  Particle  Velocity  at  N=350  -NoAR 
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Figure  4.80:  Velocity  Distributions  for  DSMC2A  Case  in  Stagnation  Region 


are  traveling  faster  as  one  would  expect.  Observing  Figures  3.4  and  3.5,  this  sample 
location  is  within  the  area  where  the  boundary  layer  and  shock  layer  merge,  so  one 
would  expect  the  particle  velocities  to  be  slower  and  for  the  distribution  to  be  fairly 
wide. 

The  differences  between  the  e  =  100%  and  NoAR  distributions  and  Bird’s  dis¬ 
tribution  are  more  noticeable  at  this  sample  location.  The  tops  of  these  distributions 
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(a)  Unweighted  EpsOO 


(b)  Weighted  EpsOO 


Unweighted  Particle  Velocity  at  N=495  -Eps50 


Weighted  Particle  Velocity  at  N=495  -Eps50 


(c)  Unweighted  Eps50 


(d)  Weighted  Eps50 


Figure  4.81:  Velocity  Distributions  for  DSMC2A  Case  Before  Step 


are  wide,  and  there  is  no  clear  point  along  the  x-axis  to  name  as  Cmp.  This  sampling 
cell  would  have  many  particles  at  varying  speeds  due  to  the  shock  layer  and  boundary 
layer  merging. 

The  sampling  location  for  Figure  4.83  is  similar  to  the  the  previous  location. 
The  sampling  cells  are  at  the  same  x  location,  but  cell  795  is  2  mm  higher  than  cell 
495.  At  this  height,  the  cell  is  not  at  the  merging  location  of  the  boundary  layer  and 
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Figure  4.82:  Velocity  Distributions  for  DSMC2A  Case  Before  Step 


shock  layer,  but  it  is  still  in  the  shock  layer,  which  means  the  particles  are  moving 
slower  than  the  bulk  velocity.  The  distributions  at  this  point  are  much  wider  than 
have  been  observed  at  the  previous  sampling  locations.  Additionally,  the  difference 
between  the  weighted  e  =  50%  distribution  and  Bird’s  distribution  is  greater  than 
previously  seen. 
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(a)  Unweighted  EpsOO 


(b)  Weighted  EpsOO 


Unweighted  Particle  Velocity  at  N=795  -Eps50 


Weighted  Particle  Velocity  at  N=795  -Eps50 


(c)  Unweighted  Eps50 


(d)  Weighted  Eps50 


Figure  4.83:  Velocity  Distributions  for  DSMC2A  Case  at  Leading  Edge  of  Step 


In  Figure  4.84,  the  weighted  and  unweighted  e  =  100%  and  NoAR  distributions 
are  similar  to  Bird’s  distribution.  The  NoAR  distributions  show  slightly  fewer  par¬ 
ticles  traveling  at  or  above  the  bulk  velocity  than  Bird.  Recalling  the  1-dimensional 
distributions,  it  should  not  be  a  surprise  that  DSMC2A  would  predict  more  particles 
traveling  at  the  bulk  velocity  compared  to  NoAR. 
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(a)  Unweighted  EpslOO 


(b)  Weighted  EpslOO 


(c)  Unweighted  NoAR 


(d)  Weighted  NoAR 


Figure  4.84:  Velocity  Distributions  for  DSMC2A  Case  at  Leading  Edge  of  Step 


Just  after  the  step  in  Figure  4.85,  all  of  the  distributions  are  very  similar  to 
Bird,  which  may  indicate  that  at  this  sampling  location  the  particles  are  near  the 
after  shock  equilibrium  state. 

Interestingly,  the  distribution  in  Figure  4.86  that  shows  the  greatest  change 
from  Bird  is  the  unweighted  e  =  100%  case. 


Ill 


(a)  Unweighted  EpsOO 


(b)  Weighted  EpsOO 


Unweighted  Particle  Velocity  at  N=1705  -Eps50 


Weighted  Particle  Velocity  at  1X1= 1705  -Eps50 


(c)  Unweighted  Eps50 


(d)  Weighted  Eps50 


Figure  4.85:  Velocity  Distributions  for  DSMC2A  Case  1mm  After  the  Step 


Farther  downstream  of  the  step,  at  sample  cell  1755,  the  particles  distributions 
are  fairly  similar.  None  of  the  distributions  are  very  smooth,  which  indicates  a  small 
sample  size.  Figure  4.56  shows  that  the  number  density  behind  the  shock  along  the  top 
of  the  step  is  low,  which  causes  the  jagged  distributions  and  causes  the  inconsistencies 
noted  in  the  contour  plots  with  Bird’s  simulations. 
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Unweighted  Particle  Velocity  at  N=1705  -EpslOO 


Weighted  Particle  Velocity  at  1X1= 1705  -EpslOO 
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(c)  Unweighted  NoAR 


(d)  Weighted  NoAR 


Figure  4.86:  Velocity  Distributions  for  DSMC2A  Case  1mm  After  the  Step 


The  weighted  e  =  100%  case  shows  the  greatest  difference  compared  to  Bird  at 
the  current  sample  location. 
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(a)  Unweighted  EpsOO 


(b)  Weighted  EpsOO 


Unweighted  Particle  Velocity  at  N=1755  -Eps50 


Weighted  Particle  Velocity  at  1X1= 1755  -Eps50 


(c)  Unweighted  Eps50 


(d)  Weighted  Eps50 


Figure  4.87:  Velocity  Distributions  for  DSMC2A  Case  Downstream  of  the  Step 


The  location  of  the  shock  layer  differed  amongst  the  simulations,  as  was  seen  in 
the  contour  plots.  Bird  and  e  =  0%  have  a  shock  layer  than  is  farther  from  the  body 
and  thicker  than  the  other  simulations.  Therefore,  Bird  and  e  =  0%  distributions  are 
very  thin  with  a  few  particles  in  the  wings.  The  unweighted  e  =  50%  distribution  is 
thin  as  well,  but  does  not  have  as  many  particles  in  the  wings,  indicating  that  the 
sample  cell  in  the  e  =  50%  is  farther  from  the  shock  layer  compared  to  Bird.  The 
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Unweighted  Particle  Velocity  at  N=1 755  -EpslOO  Weighted  Particle  Velocity  at  N=1755  -EpslOO 


(a)  Unweighted  EpslOO 


(b)  Weighted  EpslOO 


Unweighted  Particle  Velocity  at  N=1755  -NoAR 


Weighted  Particle  Velocity  at  1X1= 1755  -NoAR 


Figure  4.88:  Velocity  Distributions  for  DSMC2A  Case  Downstream  of  the  Step 


weighted  e  =  50%  has  about  the  same  number  of  particles  in  the  wings  because  these 
particles  that  are  not  at  the  bulk  velocity  have  collided,  hence  the  different  velocities. 


The  e  =  100%  and  NoAR  distributions  again  show  more  particles  near  the  bulk 
velocity,  indicating  that  the  sample  cell  is  not  in  the  shock  layer. 
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Unweighted  Particle  Velocity  at  N=3669  -EpsOO 


Weighted  Particle  Velocity  at  N=3669  -EpsOO 
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Figure  4.89: 


Velocity  Distributions  for  DSMC2A  Case  in  Shock  Layer 
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Figure  4.90:  Velocity  Distributions  for  DSMC2A  Case  in  Shock  Layer 
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Table  4.6:  Distribution  Properties  For  DSMC2A  Simulations 


Bird 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

1.9287 

0.6757 

728.18 

930.52 

795 

1.6312 

0.3718 

1004.23 

1206.81 

350 

2.4384 

1.0009 

485.69 

595.24 

3669 

7.6340 

2.3864 

2037.71 

2015.08 

1705 

1.7028 

0.5148 

1650.06 

1647.81 

1755 

1.7987 

0.3291 

1101.59 

1266.08 

EpsOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

1.9287 

0.6757 

728.18 

930.52 

795 

1.6312 

0.3718 

1004.23 

1206.81 

350 

2.4384 

1.0009 

485.69 

595.24 

3669 

7.6340 

2.3864 

2037.71 

2015.08 

1705 

1.7028 

0.5148 

1650.06 

1647.81 

1755 

1.7987 

0.3291 

1101.59 

1266.08 

Eps25 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

1.9652 

0.6749 

743.82 

924.68 

795 

1.5995 

0.3164 

1233.86 

1197.04 

350 

2.5125 

1.0360 

492.20 

596.26 

3669 

7.9278 

2.4592 

2027.63 

2021.17 

1705 

1.7581 

0.5297 

1565.48 

1658.40 

1755 

1.7277 

0.5111 

1258.97 

1263.96 

Eps50 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

2.0369 

0.7172 

697.99 

922.54 

795 

1.6470 

0.3745 

1313.31 

1187.32 

350 

2.1418 

0.8453 

514.38 

593.24 

3669 

9.6403 

2.7690 

2014.41 

2026.38 

1705 

1.8692 

0.6067 

1648.15 

1652.61 

1755 

1.5609 

0.2717 

1186.09 

1266.50 

Eps75 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

2.0761 

0.7281 

741,33 

910.85 

795 

1.5762 

0.3321 

1014.61 

1181.70 

350 

2.5055 

1.0185 

547.18 

593.34 

3669 

7.0833 

2.3225 

2035.59 

2027.47 

1705 

1.8372 

0.5897 

1734.23 

1666.23 

1755 

1.8564 

0.4565 

1358.77 

1281.05 

EpslQO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

2.0077 

0.7222 

750.70 

911.47 

795 

1.7353 

0.3966 

949.55 

1186.89 

350 

2.3224 

0.9447 

553.70 

589.91 

3669 

10.0339 

2.8572 

2029.40 

2027.33 

1705 

1.9272 

0.6783 

1530.80 

1666.63 

1755 

1.7879 

0.4319 

1258.85 

1279.68 

NoAR 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

1.8315 

0.6180 

1011.21 

916.90 

795 

1.7896 

0.5803 

980.27 

1170.35 

350 

2.1030 

0.8132 

499.81 

588.20 

3669 

7.5541 

2.4122 

2037.80 

2029.24 

1705 

1.8650 

0.6038 

1627.75 

1669.07 

1755 

1.8153 

0.5213 

1306.17 

1276.72 

EpsOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

1,9198 

0.6349 

739.12 

979.30 

795 

1,6580 

0.3800 

919.92 

1236.47 

350 

2.3359 

0.9362 

547.20 

626.53 

3669 

5.5671 

1.9010 

1995.51 

2012.54 

1705 

1.7426 

0.5000 

1661.69 

1660.13 

1755 

2.1695 

0.6175 

1092.15 

1282.78 

Eps25 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

1,9652 

0.6418 

733.16 

970.41 

795 

1,5995 

0.3577 

1222.33 

1231.20 

350 

2.5125 

0.9910 

501,24 

625.18 

3669 

7.9278 

2.0934 

2003.65 

2017.97 

1705 

1.7581 

0.5030 

1570.70 

1664.17 

1755 

1,7277 

0.6966 

1223.01 

1290.67 

Eps50 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

2.0369 

0.6681 

818.37 

969.55 

795 

1.6470 

0.4093 

1301.79 

1217.33 

350 

2.1418 

0.8041 

490.44 

623.25 

3669 

9.6403 

2.5059 

2023.47 

2023.36 

1705 

1.8692 

0.5588 

1636.00 

1664.44 

1755 

1,5609 

0.3811 

1415.57 

1295.48 

Eps75 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

2.0761 

0.6781 

773.17 

958.70 

795 

1,5762 

0.3721 

933.28 

1212.79 

350 

2.5055 

0.9395 

609.01 

621.01 

3669 

7.0833 

2.1991 

2028.38 

2028.17 

1705 

1.8372 

0.5943 

1722.11 

1682.86 

1755 

1.8564 

0.5162 

1262.01 

1303.50 

EpslOO 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

2.0077 

0.6963 

869.19 

956.80 

795 

1.7353 

0.3854 

1212.05 

1221,74 

350 

2.3224 

0.8836 

442.98 

621.19 

3669 

10.0339 

2.7731 

1982.18 

2023.66 

1705 

1.9272 

0.6393 

1518.13 

1674.38 

1755 

1,7879 

0.5120 

1376.14 

1301,92 

NoAR 

Kurtosis 

Skewness 

Cmp 

Avg  Vel 

495 

1,8315 

0.6175 

697.77 

915.52 

795 

1.7896 

0.5835 

967.21 

1172.81 

350 

2.1030 

0.8128 

492.10 

588.20 

3669 

7.5541 

2.3963 

2030.43 

2029.13 

1705 

1,8650 

0.6079 

1615.81 

1671.19 

1755 

1.8153 

0.5941 

1294.91 

1282.84 

(a)  Unweighted  °  (b)  Weighted 


In  Figure  4.6,  Bird  and  the  unweighted  e  =  0%  values  are  exactly  the  same,  as 
expected  since  the  distributions  are  colocated.  The  Cmp  value  is  the  lowest  at  sample 
cell  350  for  all  of  the  weighted  and  unweighted  simulations,  which  makes  sense  since 
the  cell  is  in  the  stagnation  region.  The  largest  Cmp  value  is  at  sample  cell  3669,  which 
is  in  the  free  stream  for  some  simulations,  but  in  the  beginning  of  the  shock  later  for 
Bird  and  e  =  0%.  Kurtosis  is  highest  for  sample  cell  3669,  as  expected  since  it  has 
the  thinnest  distribution  compared  to  the  other  cells.  The  lowest  amount  of  kurtosis 
is  located  in  cell  795,  which  was  noted  to  have  a  wide  distribution.  Cells  795  and 
1755  have  low  skewness  values,  indicating  that  the  distributions  are  more  symmetric 
compared  to  distributions  from  other  cells.  In  general,  the  weighted  distributions 
have  lower  skewness  values  compared  to  the  unweighted  distributions,  just  as  was 
seen  with  the  normal  shock  simulations. 

4-5  Hard  Sphere  Comparisons 

The  last  piece  of  this  investigation  is  determining  what  results  the  SAR  and 
NoAR  algorithms  give  when  using  the  HS  model.  Given  the  Mach  dependency  as¬ 
sociated  with  both  models,  the  question  has  been  asked  if  the  SAR  algorithm  could 
replace  the  VHS  model,  or  if  it  is  best  used  to  augment  it.  The  HS  model  was  im¬ 
plemented  for  both  the  DSMC1S  and  DSMC2A  codes,  and  the  results  have  been 
compared  to  VHS  results  and  experimental  data. 

4-5.1  1- Dimensional  Shock.  The  HS  results  were  compared  to  the  VHS 

results  through  looking  at  line  plots  and  also  comparing  the  inverse  shock  thickness 
calculations. 

The  dashed  lines  are  the  HS  results  and  the  solid  lines  are  the  VHS  results.  The 
shock  is  much  thinner  for  the  HS  model  than  the  VHS  model. 

The  offset  temperatures  can  be  see  again,  but  the  HS  results  are  farther  off  than 
the  VHS  results. 
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Figure  4.91:  Normal  Shock  Density  Line  Plot 
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Figure  4.92:  Normal  Shock  Temperature  Line  Plot 

The  magnified  figure  better  shows  how  the  HS  results  for  Bird  and  NoAR  provide 
a  consistent  temperature  as  the  VHS  results,  but  the  SAR  VHS  and  HS  results  do 
not  match. 

Again,  Figure  4.94  shows  that  the  HS  model  has  a  much  thinner  shock,  but  the 
values  before  and  after  the  shock  match. 
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Figure  4.93: 


Normal  Shock  Temperature  Line  Plot  Magnifiied 
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Figure  4.94:  Normal  Shock  U-Velocity  Line  Plot 

The  shock  in  Figure  4.95,  the  inverse  shock  thickness  is  clearly  not  consistent 
with  experimental  data,  and  SAR  and  NoAR  results  are  actually  worse  than  Bird’s 
result. 
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Reciprocal  shock  thickness  in  Ar 


Figure  4.95:  Inverse  Shock  Thickness  for  HS  Model  at  Mach  9 

4-5.2  2- Dimensional  Axisymmetric  Cylinder.  The  HS  results  provide  similar 
comparisons  to  the  Davis  experimental  data,  but  HS  tends  to  exaggerate  the  change 
in  density,  which  results  in  a  long,  pointed  profile. 


Figure  4.96:  HS  Density  Profile  at  x=. 0313m 
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In  Figure  4.96  Bird  and  e  =  0%  match  closest  to  the  experimental  data,  which 
differs  from  the  VHS  case  by  50%. 


Figure  4.97:  HS  Density  Profile  at  x=. 0462m 


In  Figure  4.97,  e  =  25%  is  the  closest  to  experimental  data,  compared  to  50% 
for  the  VHS  results. 


Figure  4.98:  HS  Density  Profile  at  x=. 0495m 


123 


Bird  and  e  =  0%  are  most  comparable  to  the  experimental  data  until  y=0.029 
m.  Afterwards  NoAR  matches  most  closely.  HS  and  VHS  results  are  very  similar  at 
the  step. 


Figure  4.99:  HS  Density  Profile  at  x=. 0509m 


Again,  the  HS  model  results  in  a  very  pointed  profile,  and  Bird  and  e  =  0% 
match  the  closest  to  Davis’  data.  The  profile  farthest  downstream  is  also  extremely 


□  Davis 

Bird 

EpsOO 

|\ 

Eps25 

I'Y 

Eps50 

Eps75 

0.04 

- 

EpslOO 

Eps200 

Eps300 

—  NoAR 

0.038 

E 

5s 


0.036 


0.034 


0.032 


0.03 


_J _ L4—+I _ I _ I _ i _ I _ I _ I _ I _ I _ I _ I _ I _ I _ i _ I _ I _ I _ I _ i _ I 

0.0001  0.00010  0.0002  0.00020 

Density,  kg/mA3 


Figure  4.100:  HS  Density  Profile  at  x=. 0561m 
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pointed  and  e  =  25%  matches  most  closely.  The  HS  model  consistently  overpre¬ 
dicts  the  density  gradient  at  the  edge  of  the  boundary  layer  compared  to  the  VHS 
models,  but  actually  underpredicts  the  density  closer  to  the  wall.  The  normal  shock 
simulations  also  showed  that  the  HS  model  overpredicts  the  density  gradient. 
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V.  Conclusions 


5. 1  Experimental  Comparison 

The  current  project  has  compared  the  new  SAR  and  NoAR  algorithms  to  three 
sets  of  experimental  data:  Alsmeyer’s  inverse  shock  thickness,  velocity  distribution 
by  Holtz  and  Muntz,  and  Davis’  2-dimensional  axisymmetric  hollow  cylinder.  The 
Alsmeyer  and  Davis  comparisons  are  based  on  the  macroscopic  density  for  the  system, 
and  certain  SAR  values  do  match  the  experimental  data  better  than  Bird.  The  inverse 
shock  thickness  plots  show  that  as  Mach  number  increases,  the  e  value  must  also 
increase  in  order  to  match  the  experimental  data.  An  empirically  derived  curve  fit  was 
used  to  relate  Mach  number  to  e.  The  comparison  to  Davis’  data  also  demonstrates 
that  SAR  can  be  used  to  change  DSMC2A  results  to  best  match  the  experimental 
data.  Due  to  the  complicated  flowfield  and  the  fact  that  the  experiment  occured  at 
only  one  Mach  number,  a  value  of  e  cannot  be  clearly  labeled  as  the  appropriate  value 
for  the  simulation,  and  a  Mach  relationship  cannot  be  identified.  The  main  point 
made  from  the  inverse  shock  thickness  and  hollow  cylinder  comparisons  is  that  the 
user  defined  e  input  allows  for  control  over  the  flowfield  in  order  to  best  match  the 
experimental  data.  The  potential  exists  after  further  research  for  SAR  to  be  used  as 
method  of  producing  more  accurate  results  for  a  multitude  of  DSMC  simulations. 

The  velocity  distributions  allow  for  a  deeper  investigation  into  the  effect  of 
SAR  and  NoAR  on  the  flowfield.  Holtz  and  Muntz  performed  an  experiment  at  Mach 
7.18  to  find  the  velocity  distributions  of  the  particles  through  a  normal  shock  layer. 
Comparing  these  results  to  Bird,  SAR,  and  NoAR  shows  that  in  equilibrium,  the 
distributions  are  very  nearly  the  same.  The  idea  that  in  equilibrium  the  algorithms 
produce  similar  results  is  confirmed  when  looking  at  the  macroscopic  results.  For  the 
1 -dimensional  normal  shock,  the  line  plots  show  that  before  and  after  the  shock  the 
values  are  the  same,  and  only  through  the  shock  do  the  flowfield  properties  change 
with  the  exception  of  temperature.  There  is  a  problem  with  the  way  temperature  is 
sampled  which  results  in  a  higher  value  compared  to  Bird  for  SAR,  but  not  NoAR. 
The  DSMC2A  code  samples  temperature  the  same  way,  but  does  not  suffer  from 


126 


the  same  problem.  Further  research  into  this  matter  is  required.  For  DSMC2A, 
the  percent  difference  plots  show  very  little  difference  in  flowfield  parameters  in  the 
freestream.  The  differences  were  noticed  in  the  shock  layer,  boundary  layer,  and 
stagnation  region.  These  are  areas  where  the  gradients,  and  thus  nonequilibrium 
are  high.  Velocity  distributions  were  taken  at  three  locations  within  the  shock  and 
compared  to  Holtz  and  Muntz.  The  plots  definitively  show  that  the  SAR  algorithm 
agrees  with  the  experimental  velocity  distributions  the  best.  The  SAR  algorithm 
differs  from  Bird’s  code  in  two  respects:  the  collision  rate,  and  flowfield  sampling. 
The  change  in  collision  rate  does  alter  the  flowfield,  but  not  significantly.  The  real 
change  is  due  to  the  flowfield  sampling.  SAR  and  NoAR  only  sample  particles  that 
have  collided  during  the  current  time  step.  The  partial  weighting  associated  with 
SAR  allows  particles  to  collide  that  normally  would  be  rejected.  While  these  partial 
collisions  are  weighted  less  than  one,  their  inclusion  changes  the  velocity  profiles  and 
allows  for  more  accurate  results.  Bird’s  resulting  velocity  profiles  show  a  significant 
number  of  particles  at  the  bulk  velocity  deep  in  the  shock  layer,  which  according 
to  the  experimental  data  is  not  reasonable.  The  collision  rate  does  not  allow  for 
local  equilibration,  which  is  a  key  assumption  of  DSMC,  and  the  reason  the  velocity 
distributions  show  a  spike  at  the  bulk  velocity  within  the  shock.  NoAR  tended  to 
produce  velocity  distributions  that  underestimate  the  number  of  particles  traveling 
slower  than  the  bulk  velocity  within  the  shock  layer,  as  did  e  =  100%. 

5. 2  Further  Results 

There  is  also  important  information  to  be  gleaned  from  the  other  data  not 
compared  to  experimental  data.  Velocity  distributions  were  created  for  Mach  numbers 
1,  3,  6,  and  9.  The  lower  Mach  numbers  of  1  and  3  show  that  the  DSMC  simulation 
is  able  to  maintain  the  local  equilibrium  distribution  within  the  cell.  The  gradient 
within  the  shock  is  shallow  enough  that  DSMC  is  accurate,  which  is  why  through 
Mach  3,  the  inverse  shock  thickness  values  are  accurate.  At  Mach  numbers  6  and  9, 
the  velocity  distributions  within  the  shock  are  bimodal,  indicating  that  the  particles 
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are  not  in  equilibrium  with  each  other.  The  failure  to  maintain  local  equilibrium 
causes  the  simulation  to  take  longer  to  equilibrate  through  the  shock,  which  results 
in  a  thicker  shock  compared  to  experimental  data.  One  situation  that  can  cause  a 
simulation  to  not  maintain  local  equilibrium  is  a  cell  that  is  larger  than  the  mean  free 
path.  In  a  shock  layer,  the  density  increases,  which  would  in  turn  decrease  the  mean 
free  path.  The  mean  free  path  in  the  freestream  is  an  order  of  magnitude  larger  than 
the  cell  size,  but  to  rule  out  the  possibility,  a  simulation  with  the  cell  size  reduced  by 
half  was  run.  The  resulting  velocity  distribution  in  the  shock  showed  little  appreciable 
difference  from  the  original  velocity  distribution. 

Running  the  DSMC2A  code  with  a  grid  that  was  a  quarter  the  size  of  the  original 
grid  resulted  in  requiring  a  smaller  e  value  to  match  results.  The  smaller  grid  allowed 
for  the  particles  to  maintain  local  equilibrium  better,  which  means  a  smaller  e  value 
can  be  used,  just  as  a  smaller  e  value  can  match  the  inverse  shock  thickness  for  lower 
Mach  numbers.  The  required  e  input  value  is  proportional  to  Mach  number  because 
the  higher  Mach  numbers  create  more  nonequilibrium  within  the  shock  layer,  which 
needs  to  be  overcome  by  a  larger  e  to  allow  for  proper  equilibration. 

The  VHS  algorithm  is  already  known  to  produce  better  results,  and  Figure  4.95 
shows  that  even  with  SAR,  the  HS  results  do  not  compare  well  to  experimental  data. 
The  HS  density  plots  in  Figures  4.96  through  4.100  show  an  exaggerated  change  in 
density  compared  to  the  VHS  density  profiles,  and  in  general  the  SAR  and  NoAR 
results  are  actually  farther  from  experimental  data.  Therefore,  the  HS  model  cannot 
be  used  in  conjunction  with  SAR  to  produce  accurate  results. 

5. 3  Future  Work 

For  SAR  to  be  used  in  simulations  a  few  things  need  to  first  occur:  SAR  needs 
to  be  compared  to  more  complex  simulations  and  geometry,  and  SAR  needs  to  be 
compared  to  experiments  having  higher  Mach  numbers  than  the  Alsmeyer  data  con¬ 
tains.  In  order  for  SAR  to  be  used  on  cases  with  more  complexity,  SAR  will  need  to 
be  applied  to  a  DSMC  code  that  allows  for  geometries  such  as  a  blunt  body,  which  has 
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been  used  in  experiments.  Bird  has  not  made  his  more  recent  source  codes  available 
to  the  public,  so  a  DSMC  code  would  need  to  be  developed  specifically  with  the  SAR 
modifications,  or  another  open  source  DSMC  code  could  be  used  for  modification. 
The  results  of  the  SAR  simulations  can  be  compared  to  experimental  data.  Addition¬ 
ally,  in  order  to  understand  the  Mach  dependency,  further  experimental  comparisons 
that  vary  Mach  number  throughout  the  hypersonic  regime  will  need  to  be  completed. 
Additionally,  the  DSMC2A  code  was  compared  to  experimental  data  at  one  Mach 
number,  data  that  varies  Mach  number  should  be  used  to  best  compare  DSMC  to 
SAR  and  NoAR  implemented  into  a  2-dimensional  algorithm. 
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